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1.  Introduction 


The  process  of  obtaining  images  of  the  reflectivity  or  density  of  target 
areas  that  are  rotating  and  translating  with  respect  to  a  sensor  such  as  a 
mono-static  or  bistatic  radar,  a  sonar,  or  an  x-ray  CAT  scanner,  has  been 
studied  for  the  past  40  years  1*2,3, 4  Ausherman  etal.5  have  written  an 
excellent  review  of  the  work  done  in  this  area.  Target  areas  that  rotate 
and  translate  relative  to  a  radar  include,  for  example,  planet  surfaces 
observed  from  a  satellite,  ground  terrain  observed  from  air- borne 
platforms,  sub-surface  objects  and  voids  observed  from  moving  vehicles, 
and  people  scanned  by  a  rotating  x-ray  system.  Although  the  processing 
described  is  applicable  to  other  systems,  this  article  shall  treat  the  topic 
from  the  point  of  view  of  a  SAR  (Synthetic  Aperture  Radar).  As  the 
aspect  angle  between  the  sensor  and  a  target  changes  with  time,  the  sensor 

collects  a  sequence  of  signal  records.  After  collecting  data  for  Ts  seconds, 

the  aspect  angle  has  changed  by  6S  degrees.  These  received  signals  are 
then  coherently  processed  to  obtain  the  reflectivity  profile  of  the  target 
area.  Down-range  resolution  into  range-bins  is  obtained  primarily  by  the 
bandwidth  of  the  sensor.  Cross-range  resolution  is  obtained  primarily  by 
coherently  processing  the  received  signals  such  that  a  very  wide  aperture 

is  simulated;  an  aperture  that  is  6S  degrees  wide. 

Implementation  and  study  of  this  image  formation  processing  topic  has 
been  limited  in  two  ways.  First,  the  image  formation  processing  has 
assumed  that  targets  are  isotropic  point  scatterers.  Many  targets,  however, 
are  anisotropic  resonant  scatterers.  Treatment  of  this  resonant  case 
becomes  important  when  the  sensor  spectrum  covers  the  Rayleigh, 
resonant,  and  optical  regions  of  a  family  of  targets.  For  example, 
discrimination  between  scatterers  can  be  based  on  the  unique  signature  of 
each  target.  But,  in  order  for  the  discrimination  to  work  in  the  context  of 
microwave  reflectivity  imaging,  the  information  in  the  signature  must  be 
preserved  during  the  image  formation  process. 


*  William  M.  Brown  and  Ronald  J.  Fredricks,  "Range-Doppler  Imaging  with  Motion  through  Resolution  Cells", 
IEEE  Trans.  Aerosp.  Electron.  Syst.,  AES-5  No.l,  (Jan  1969),  98-102 

2 J.L .  Bauck  and  W.K.  Jenkins,  'Tomographic  Processing  Of  Spotlight-Mode  Synthetic  Aperture  Radar  Signals 
With  Compensation  For  Wavefront  Curvature",  IEEE  Interna.  Conf.  on  Acoust.  Speech  and  Sig.  Proc.  (Apr. 11-14 
1988 ) 

3 Edwin  D.  Banta,  "Limitations  on  SAR  Image  Area  Due  to  Motion  Through  Resolution  Cells",  Correspondence, 
IEEE  Trans.  Aerosp.  Electron.  Syst.,  AES-22,  No.  6  (Nov  1986),  799-803 

^D.  Mensa,  and  G.  Heidbreder,  "Bistatic  Synthetic-Aperture  Radar  Imaging  of  Rotating  Objects",  IEEE  Trans. 
Aerosp.  Electron.  Syst.,  AES-18,  No.  4  (July  1982),  423-431 

5 Dale  A.  Ausherman,  Adam  Kozma,  Jack  L.  Walker,  Harrison  M.  Jones,  and  Enrico  C.  Poggio,  "Developments  in 
Radar  Imaging”  IEEE  Trans.  Aerosp.  Electron.  Syst.  AES-20,  No. 5  (July  1984),  363-400 
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The  second  limitation  in  previous  work  is  that  relatively  narrow 
bandwidth  systems  have  been  assumed.  For  example,  the  compressed 
pulse  width  of  the  sensor  is  assumed  to  be  at  least  several,  and  usually 
many  cycles  of  an  RF  carrier  frequency.  So  a  single  range-bin  is  derived 
from  several  cycles  of  the  carrier.  Newer  UWB  (Ultra  Wide  Bandwidth) 
sensors,  however,  have  made  it  possible  to  make  the  image  range-bin  size 
roughly  1/2  the  wavelength  of  the  highest  frequency  in  the  sensors’ 
spectrum.  The  bin  size  is  much  smaller  (1/10  or  less)  than  the  wavelength 
of  the  lowest  frequency  in  the  spectrum.  This  wide  bandwidth 
exacerbates  range  walk  and  wavefront  curvature  errors  to  the  point  where 
conventional  FFT  based  processing  must  be  restricted  to  very  small 
patches  within  an  image  area. 

These  two  limitations  are  not  independent.  It  is  bandwidth  that  allows 
target  discrimination  based  on  signature  analysis  and  it  is  bandwidth  that 
makes  image  formation  more  difficult.  It  is  the  purpose  of  this  paper  to 
examine  image  formation  processing  that  preserves  resonant  target 
signatures  and  to  present  an  efficient  method  of  solving  the  microwave 
reflectivity  imaging  problem  for  UWB  signals  and  resonant  targets. 


6 


2.  Background 

2.1.  Terminology 

SAR  systems  depend  upon  collecting  data  coherently  along  a  path.  This 
path  is  referred  to  as  the  "synthetic  aperture."  Figure  1  shows  a  sketch  of 
the  scenario  projected  onto  a  plane.  An  aircraft  flies  at  some  elevation  h 

above  the  ground,  over  a  distance  Ls  to  form  the  synthetic  aperture.  The 
image  area  grid  is  referenced  to  the  center  of  the  aperture.  Range  is 
marked  off  by  the  parameter  i.  The  parameter  k  marks  off  azimuthal 
(bearing)  lines.  The  sample  points  in  the  aperture  are  marked  off  by  the 
parameter  j.  The  distance  between  the  /h  point  in  the  aperture  and  the 

fl | 

(i,  k)  position  in  the  image  area  is  denoted  by  di  j  k .  Although  nearly  all 

operational  SAR  systems  use  a  rectangular  grid,  there  are  advantages  to 
the  polar  grid  that  will  be  discussed  later  in  this  report.  Both  the 
azimuthal  lines  and  the  range  lines  are  referenced  to  the  center  of  the 
aperture.^1  Table  1  summarizes  the  nomenclature  to  be  used  in  the  rest  of 
the  report. 


i.K)  -bt-2 


Figure  1.  Basic  Flight  Path  and  Image  Area  Geometry 


®This  grid  is  convenient  for  post  processing.  An  X-Y  grid  can  also  be  computed  using  the  techniques  described. 
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Table  1.  Analysis  Nomenclature 


4 

[Length  (meters)  of  the  synthetic  aperture. 

Sj(t) 

Denotes  a  received  signal  amplitude  (volts)  as  a  function  of  time  (seconds 
from  the  leading  edge  of  the  transmitted  pulse),  at  the  jth  position  in  the 
aperture.  The  A/D  outputs  represent  this  signal. 

Over  type  hat 
e.g.  s(t) 

Denotes  that  the  function  is  an  estimate.  The  example  denotes  an  estimate  for 
a  received  signal. 

Subscript  of 
i,kj 

Particularizes  the  function  to  be  at  particular  geometry's,  like  the  position 
in  the  synthetic  aperture,  or  the  bearing  line,  or  the  (i,/:)^  position  in  the 

image  area. 

Denotes  the  distance  (meters)  between  the  j^1  position  of  the  radar,  and  the 
(i,k)th  position  in  the  area  to  be  imaged. 

El 

Denotes  the  round  trip  time  (seconds)  for  the  radio  energy  to  travel  from  the 
jth  position  in  the  synthetic  aperture,  to  the  (i,k)th  position  in  the  image  area, 
and  back. 

L  c 

The  speed  of  light. 

I 

Time  shift  between  range  bin  i  and  range  bin  i+1  at  the  center  of  the  aperture 
(where  ./=()). 

| 

Pulse-Repetition-Frequency  of  the  radar  in  Hertz. 

V- 

Relative  Bandwidth  u  -  (Fhi  -Fh)/  F0  and F0  =  (Fw  +  Fh ) / 2 

UWB 

Ultra-Wide  Bandwidth,  where  jl~l 

X 

Wavelength  in  meters. 

2.2.  Approaches  to  Focusing 

SAR  processing,  or  SAR  focusing,  is  sometimes  referred  to  as  a  Doppler 
based  process.  Why  is  this  "Doppler"  paradigm  used?  Referring  to 
Figure  1,  assume  that  the  radar  platform  is  moving  at  velocity  v  as  it 
collects  data  along  the  aperture.  Therefore,  the  data  collection  points  are 
spaced  equally  in  time,  occurring  at  a  rate  governed  v  and  the  PRF.  (This 
sampling  across  the  aperture  is  sometimes  referred  to  as  "slow-time".)  If 

the  radar  is  operating  at  a  frequency  f(),  then  the  echo  from  a  target  in  the 
image  area  will  have  a  Doppler  shift  profile  as  the  platform  moves 
through  the  aperture.  It  will  have  an  upward  shift  while  the  platform 
approaches,  the  shift  will  drop  to  zero  as  the  platform  moves  to  a  position 
broadside  to  the  target,  and  the  shift  will  be  downward  as  the  platform 
recedes.  Every  target  position  will  have  a  unique  Doppler  profile.  Since 
every  position  in  the  image  has  a  distinct  Doppler  profile,  an  image  can  be 
formed  by  simply  assigning  to  each  pixel  a  filter  matched  to  the  Doppler 
profile  expected  for  a  target  at  the  location  represented  by  that  pixel. 
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The  classic  "Polar  Formatting"^  approach  to  computing  a  SAR  image 
simply  adjusts  (time  shifts)  the  data  at  each  aperture  point  so  that  the  new 
data-set  appears  as  if  the  platform  had  moved  in  a  short  circular  path  - 
with  the  circle  centered  at  the  center  of  the  image  area.  Once  this 
formardng  is  done,  a  target  at  the  center  point  has  the  unique  Doppler 
profile  of  zero  over  the  entire  aperture.  Other  points  have  other  unique 
Doppler  profiles.  A  Fourier  transform  is  typically  used  to  recover  the 
image  from  the  Doppler  profiles.  The  Fourier  approach  works  well  as 
long  as  the  circular  path  is  sufficiently  short,  the  image  area  sufficiently 
small,  the  signal  bandwidth  sufficiently  small,  and  the  distance  from  the 
radar  to  the  image  center  sufficiently  long.  This  paper  addresses  the  case 
where  none  of  these  restrictions  apply. 

Although  the  Doppler  paradigm  has  helped  countless  people  to  visualize 
SAR  focusing,  it  is  not  necessary,  nor  always  helpful,  in  solving  the  SAR 
focusing  problem.  For  example,  Doppler  shift,  which  is  defined  as  2v/X, 
works  fine  when  X  varies  by  a  few  percent,  but  it  becomes  a  stumbling 
block  in  the  UWB  case  where  X  can  vary  by  10  or  100  to  1.  Doppler  has 
proven  to  be  a  very  convenient  narrow  band  concept.  But  its  convenience 
breaks  down  at  wide  relative  bandwidths.  The  focusing  problem  can  also, 
however,  be  looked  at  as  a  stationary  array  of  N  antennas  whose  outputs 
are  digitally  stored  and  combined  in  a  computer  to  form  beams.  It  is  the 
author's  view  that  the  focusing  problem  is  easier  to  visualize  and  solve 
using  this  stationary  array  paradigm  for  the  case  where  neither  the 
geometry  nor  the  bandwidth  are  restricted. 

Using  the  stationary  array  paradigm  one  can  see  that  equation  1  coherently 
focuses  the  SAR  data  by  summing  across  the  array.  It  may  be  helpful  to 
point  out  that  both  approaches  can  be  seen  to  be  identical  at  the  center 
point  of  the  "Polar  Format"  patch.  Notice  that  regardless  of  carrier 
frequency,  the  center  point  of  the  polar  formatted  data  is  always  analyzed 
by  the  DC  term  of  the  Fourier  transform.  And  the  DC  term  of  a  Fourier 
transform  is  simply  a  summation  of  the  data  points.  The  summation 
shown  in  equation  1  is  identical  to  finding  the  DC  term  of  the  polar 
formatted  SAR  data.  But  instead  of  formatting  once  and  then  finding 
many  "Doppler”  profiles  via  an  FFT,  equation  1  formats  and  "sums”  the 
data  many  times;  each  formatting  makes  a  different  pixel  the  center  and 
then  the  DC  term  summation  is  calculated  to  get  the  value  for  that  pixel. 
The  result  is  that  optimally  focused  beams  are  formed  -  beams  taking  full 
advantage  of  both  the  entire  aperture  and  the  entire  signal  bandwidth.  The 
focusing  is  truly  frequency-independent. 


7j.  L.  Walker,  "Range-Doppler  imaging  of  rotating  objects,"  IEEE  Trans.  Aerosp.  Electron.  Syst.,  AES-16  (Jan. 
19S0),  23-52 
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How  does  the  beam  width  compare  with  "conventional"  (far-field  narrow 
band)  antenna  theory?  The  rule-of-thumb  half-power  beamwidth  of  a  line 
array  is  approximately  2JL  where  L  is  the  length  of  the  array.  The  beams 
formed  by  the  processing  described  above  follow  this  rule  and  have  a 
width  proportional  to  A.  Low  frequencies  have  wide  beams  and  high 
frequencies  have  narrow  beams.  An  interesting  aspect  of  this  fact  is  how 
it  manifests  itself  in  the  time-domain  as  a  source  moves  through  a  beam. 
The  impulse  response  at  the  center  of  the  beam  is  a  nice  narrow  pulse. 
But  as  one  moves  away  from  the  center  of  the  beam,  the  impulse  response 
gets  broader  and  broader.  This  time-domain  broadening  occurs  because 
more  and  more  high  frequency  energy  is  lost  as  the  source  moves  out  of 
the  narrowing  high-frequency  beam. 

As  soon  as  one  switches  from  the  Doppler  to  the  stationary  array 
approach,  it  may  be  tempting  to  also  switch  to  thinking  in  "phased  array" 
terms.  To  do  so  is  a  mistake  when  bandwidth  and/or  geometry  are  not 
restricted.  Whenever  the  geometry  is  not  restricted  to  the  far-field  of  an 
aperture,  plane-wave  simplifications  break  down.  This  break  down 
invalidates  simple  phase  steering.  Since  a  Fourier  transform  forms  beams 
by  simple  phase  steering,  Fourier  techniques  becomes  less  and  less  useful 
as  targets  move  into  the  near-field.  In  an  ultra  wide  bandwidth  system, 
phase  becomes  meaningless  with  regard  to  defining  element  positions  or 
the  beam  forming  network.  To  speak  of  shifting  one  element  180  degrees 
with  respect  to  another  element  implies  a  fixed  A.  If,  for  example,  A 
changed  by  2  to  1,  then  a  delay-line  that  provided  180  degrees  at  one 
frequency  would  provide  360  degrees  at  the  other.  Yet  delay-lines  are 
precisely  the  element  needed  to  build  a  wide  bandwidth  "phased  array" 
antenna.  When  A  varies  several  octaves,  the  best  parameters  to  use  in  the 
equations  defining  the  antenna  is  the  time  and  distance  -  time-shift  of  the 
delay  lines  that  form  the  beam-forming  network,  and  distance  between 
antenna  elements.  So  it  is  best  if  one  switches  to  thinking  in  "timed  array" 
terms.  This  time-based  framework  results  in  derivations  that  are  frequency 
and  geometry  independent. 

2.3.  Null  Steering  and  Pattern  Forming 

Generally,  an  antenna  designer  would  say  that  the  most  critically 
important  aspect  of  making  desirable  antenna  patterns  is  forming  nulls.  A 
simple  classic  case  is  spacing  2  elements  at  90°  and  phasing  them  by  90° 
to  form  an  endfire  beam  in  one  direction  and  a  null  in  the  opposite 
direction.  In  the  simplest  case  (using  element  weighting  of  +1  and  -1),  to 
form  a  UWB  null  one  must  time-steer  the  array  to  where  the  null  should 
be,  and  then  invert  half  the  elements  prior  to  summing.  Of  course,  using 
other  weighting  factors  allows  more  freedom.  An  impulse  signal  8(t) 
coming  from  a  direction  other  than  the  null  would  produce  in  the  receiver 
some  array-induced  waveform.  Ignoring  bandwidth  effects  from  each 
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element,  that  waveform  would  be  a  function  of  the  +  and  -  weighting  and 
element  spacing.  And  by  definition,  that  waveform  is  the  antenna-array 
impulse  response  for  that  beam  angle.  Matched  filtering  to  that  waveform 
forms  a  mainlobe  in  that  direction. 

Grating  lobes  are  customarily  defined  as  lobes  whose  gain  is  equal  to  the 
mainbeam.  It  is  interesting  to  note  that  this  definition  leads  to  confusion 
in  the  UWB  case.  Are  there  grating  lobes?  Well,  that  depends.  If  one 
looks  at  the  response  to  the  UWB  signal,  like  an  impulse,  then  the  answer 
is  no.  The  peak  response  on  the  mainlobe  is  higher  than  the  response  at 
any  other  angle.  So  the  definition  fails.  On  the  other  hand,  however,  the 
antenna  is  a  linear  system.  It  behaves  just  like  an  identical  array  operating 
at  a  single  frequency.  So,  if  one  looks  at  the  response  to  a  CW  signal,  the 
answer  can  be  yes.  If  the  elements  are  physically  spaced  greater  than  A/2 
apart  at  the  highest  frequency  component  in  the  UWB  waveform,  then 
yes,  there  certainly  will  be  a  grating  lobe  at  that  frequency  component. 
All  the  insight  gained  from  CW  antenna  analysis  holds  and  remains 
useful.  One  must,  however,  be  careful  when  applying  terms  like  grating- 
lobes  that  presuppose  a  narrow  band  signal.  In  the  case  of  SAR,  the 
element  spacing  (Velocity/PR  F)  needed  will  be  a  function  of  what 
sidelobes  are  permissible  at:  the  various  frequency  components  in  the 
UWB  waveform. 

2.4.  Target  Resonance  Effects 

Historically,  the  relative  bandwidth  of  radars  has  been  sufficiently  small 
that  a  target's  echo  is  adequately  modeled  by  a  single  number,  [sigma],  the 
Radar  Cross  Section  (RCS);  usually  given  in  square  meters.  When  pi  >.5, 

however,  a  single  number  may  no  longer  adequate  --  a  is  a  function  of 
frequency.  For  example,  Figure  2  is  a  plot  of  the  RCS  of  a  sphere.  In 
addition  to  the  magnitude  characteristic  plotted,  there  is  also  a  phase 
characteristic.  These  two  frequency-domain  characteristics  can  also  be 
represented  in  the  time  domain  by  a  ringing  or  resonant  response.  In 
either  case  -  time  domain  or  frequency  domain  -  the  plots  can  be  referred 
to  as  the  impulse  response  of  the  target. 
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2 Ttfr  /  c 

Figure  2.  Radar  Cross  Section  of  a  Sphere 

Since  historical  SAR  systems  have  not  been  interested  in  target  ringing, 
no  attention  was  paid  on  how  to  preserve  the  ringing  information.  This 
paper  describes  and  analyzes  a  procedure  to  focus  a  large  array,  over  ultra¬ 
wide  bandwidths,  in  such  a  way  as  to  preserve  the  resonant  response  of 
targets;  even  when  they  are  in  the  near-field. 
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3.  Fundamental  Time  Based,  Physical  Array  Approach 

Consider  an  aperture  looking  at  completely  empty  space  except  for  an 
isotropic  scatterer  at  position  (i,k).  Also  assume  that  an  ideal  impulse  <$(/) 
is  broadcast  It  is  desirable  to  take  advantage  of  all  the  echo  energy  across 
the  aperture.  To  do  that,  the  energy  from  all  collection  points  along  the 
synthetic  aperture  is  summed.  However,  this  summation  must  be  done 
such  that  the  energy  adds  coherently  at  all  frequencies.  Frequency 
independent  adding  is  accomplished  as, 

fi,k  (0  =  X sj WiXj  +  *>  for  t  <  0  (1) 

; 

Here,  fy^t)  would  be  the  impulse  response  of  the  target  located  at  position 
Note  that  the  Tikj  term  time  shifts  the  received  signals  sy  such  that 
the  target  impulse  response  starts  at  t= 0  —  at  all  points  in  the  aperture. 

The  above  discussion  assumes  that  a  target's  impulse  response  is 
sufficiently  similar  at  all  points  along  the  aperture  to  consider  equation  1 
true.  The  response  of  a  vertical  dipole,  for  example,  does  not  change 
depending  on  where  the  radar  is  positioned  in  the  aperture.  It  is, 
therefore,  an  isotropic  target. 

Suppose  we  relax  the  isotropic  requirement,  and  suppose  that  a  complex 
scatterer  is  at  position  (i,k)  in  the  image  area.  A  horizontal  dipole,  for 
example,  is  an  anisotropic  target  Advantage  could  be  taken  of  the 
anisotropic  behavior  to  extend  the  detection  and  target  recognition 
performance  of  the  radar  by  rewriting  equation  1  as 


fi,k (0  =  X kj  ®  SjWiXj  +  o  for  I  <  0 ,  (2) 

/' 

where  the  convolution  step  represents  filtering.  In  this  case,  a  set  of  filters 
Xj( CO)  would  be  needed,  each  one  matched  to  the  target  response  at  the 
bearing  of  the  /h  position  in  the  aperture.  In  order  to  simplify  the  rest  of 
this  report,  equation  1  will  be  used  as  the  fundamental  equation. 
Nonetheless,  one  must  recognize  that  real  targets  are  complex  anisotropic 
polarimetric  scatterers.  This  fact  should  not  be  ignored  for  large  arrays. 
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4.  Short  Impulse  Response  Approximation 

A  typical  2-D  SAR  image  is  simply  the  echo  magnitude  mapped  to 
intensity.  Equation  1  expands  the  typical  2-D  image  into  a  3-D  image 
with  the  target  ringing  along  the  third  dimension.  Given  the  heavy 
computation  load  of  typical  SAR  processing,  if  it  is  required  that  an /(f)  be 
computed  for  every  pixel  instead  of  a  single  value,  then  a  massive 
computer  would  be  needed.  The  impact  of  this  computational  load  is  even 
worse  when  one  considers  that  the  mass  of  resulting  data  must  be  analyzed 
in  the  target  detection/identification  phase.  This  section  defines  an 
approximation  which  reduces  the  problem  back  to  a  2-D  case  and  presents 
an  error  analysis  of  the  approximation. 

4.1.  Theory 

The  problem  is  that  a  third  dimension  has  been  added  to  measure  target 
ringing.  Note,  however,  that  if  there  was  only  one  aperture  position,  then 
ringing  of  the  target  would  just  appear  in  pixels  behind  the  target.  Even 
with  an  aperture  of  many  positions,  ringing  will  appear  behind  the  target. 
But  since  the  geometry  is  not  constrained,  and  near-field  operation  is 
presumed,  the  antenna  beam  defocuses  behind  the  target.  Thus  (1)  will 
perfectly  focus  ringing  of  unbounded  duration.  But  for  practical  purposes, 
the  target  only  rings  for  a  finite  duration,  say  M  range-bins.  The  question 
is,  given  that  the  ringing  is  finite  in  duration,  can  the  ring  information  be 
retained  in  a  2-D  image.  In  other  words,  suppose  that  instead  of 
calculating  a  time  series  /*(f)  for  each  pixel  in  the  image,  only  one  value 
is  calculated,  for  example:  /^(t)  where  x  is  fixed  for  the  image.  If  the 
target  identification  problem  can  be  solved  with  this  2D  image,  then  a 
great  reduction  in  computational  load  is  obtained.  The  aim  of  this  section 
is  to  describe  and  quantify  the  error  bounds  when  a  2-D  image  fi  k{ x)  is 
used. 

First,  define  a  to  be  the  round-trip  time  it  takes  the  radar  pulse  to  traverse 
one  range  bin  on  the  grid;  that  is, 

a  =  ti+i,k,i~fi,/c,j:Vk'j  =  °  (3> 

Since  the  grid  for  the  image  has  been  defined  to  be  referenced  to  the 
center  of  the  aperture  (the  ./=()  position),  the  i  parameter  can  be  thought  of 
as  a  quantized  time  t  parameter.  While  t  is  in  units  of  seconds,  i  is  in  units 
of  range  bins,  and  they  are  related  by  a  —  one  range  bin  equals  a  seconds. 
Now  we  can  write 


Sj C T^j  +ma)  =  Sj (Ti+mJcJ )  for 


[case  1:  jf'  =  0,Vm,Vi,  Vfcl 
{case  2:  m  —  0,V/,V7,V£j 
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An  approximation  is  made  to  equation  4,  generalizing  it  to  include  all 
aperture  points  over  a  limited  range  of  m  to  arrive  at 

Sj (Ti  k  j  +  ma)  »  Sj (Ti+mk  j )  for  Vi,  V£,  V/, m  =  -M- -  •  M  (5) 

The  following  may  be  said  of  the  approximation  in  Equation  5. 

1 .  It  is  perfect  at  the  center  of  the  aperture  regardless  of  m. 

2.  It  is  perfect  at  m= 0  regardless  of  j; 

3.  gets  worse  as  m  deviates  further  from  zero; 

4.  and  is  worst  at  the  end-points  of  the  aperture. 

A  solution  is  desired  to  (1)  in  the  form  offi  k(x)  where  x  is  a  constant.  In  a 
practical  system  x  will  be  mapped  to  discrete  range  bins,  so  let 


t  =  na.  (6) 

Substituting  equations  5  and  6  into  1  we  define  a  2-D  "image  "/(i,£)  as 

/w<°  >  =  l¥W 

i 

If  one  thinks  of  a  ringing  target,  then  equation  7  can  be  described  as 
allowing  the  point  of  perfect  focus  to  be  adjusted  to  any  depth  n  in  the 
ring.  For  example,  if  n= 0,  then  the  perfect  focus  point  would  be  at  the 
leading  edge  (the  first  sample)  of  the  ringing  response.  If  n  *  0,  say  n=3, 
then  the  third  sample  of  the  ringing  response  would  be  perfectly  focused. 

Next  consider  the  indexing.  The  indexing  is  performed  such  that  f(i,k) 
always  represents  the  leading  edge  of  the  response  from  a  target  located  at 
(i,k)  regardless  of  whether  it  is  perfectly  focused  or  not.  Once  f(i,k)  is 
found,  we  now  wish  to  find  the  discrete  samples  of  the  target  ringing.  The 
samples  will  be  counted  as  m=Q  for  the  first  sample,  m=  1  for  the  second 
and  so  on.  These  samples  are  obtained  by  simply  incrementing  the  i 
index.  The  value  is  just  f(i+m,k );  which  follows  from  equations  5 
and  7  as 

ft,k (ma)  ~  fi-n+mik(na)  =  f(i+  m, k)  or, 

Zs;<Tf*,/+™a)  ~fi-n+m,k(na)  =  f(i  +  tn,k)  (8) 

i 

Clearly  equation  8  is  identical  to  equation  1  (i.e.  perfect)  when  m=n.  So 
equation  8  gives  perfect  focus  at  m=n..  The  name  "short  impulse  response 
approximation"  is  used  because  the  approximation  only  needs  to  remain 
accurate  over  the  duration  of  a  target's  impulse  response. 


15 


4.2.  Examples 

To  illustrate  the  use  of  equation  8,  we  follow  these  steps: 

1.  Fix  n; 

2.  Place  a  target  (for  the  purposes  of  the  illustration)  at  say  i=: 237 
in  range  on  the  k"1  bearing. 

3.  Use  equation  7  to  calculate  f(i,k)  for  all  (i,k); 

The  cases  of  interest  are  where  n= 0  and  where  n^O.  We  will  consider 
each  separately. 

4.21.  Case  where  n=0 

Use  equation  8  to  find  the  ringing  response  of  the  target.  The  response  for 
the  first  three  points  of  the  focused  impulse  response  is: 

fi2>7 ,k  (0)  = /237-n+m,jt(0)  =  /237jA:(0)  =  /(237/fc)  (case  for  m  =  0} 
/237^(a)  -/237_n+m^(0)  =  /238/fc(0)  =  /(238^)  {case  for  m  =  l) 
fi37,k  (2a)  =  /237-H+m,fc(0)  =  /239jfc(0) =  /(239,k)  {case  for  m  =  2} 

Since  n=0,  the  amplitude  of  the  leading  edge  (m= 0)  of  the  impulse 
response  of  that  target  is  perfectly  focused.  As  m  increases,  the 
approximation  gets  worse.  So  the  approximation  is  useful  as  long  as  the 
target  resonance  dies  before  the  approximation  gets  too  bad. 

4.2.2.  Case  where  n  A  0 

In  this  case,  perfect  focus  is  na  seconds  past  the  leading  edge  of  the  target 
impulse  response.  Suppose,  for  this  example,  n-2.  The  first  4  data  points 
for  the  impulse  response  is: 

/237,*(0)  =5  /237-n+m,Jt(««)  =  /235,fc(2a,)  =  /(237,fc)  {case  for  m  =  0} 
/237,fc(«)  af237-n+mAna)  =  W  (2o:>  = /(238' ^  <case  for  m  =  1} 
fi37,k (2a)  =  /237-M+BI/jfc(««)  =  /237(fc(2a)  =  /(239,fc)  {case for m- 2} 
fi37,kOa)  =  /237-n+m/fc(»«)  =  /23s,*  <2a)  =  /( 240, k)  {case for  m  =  3} 


Note  that  the  leading  edge  is  not  perfectly  focused,  as  it  was  when  n=0. 
Incrementing.  The  important  point  here  is  that  one  can  choose  n  #  0  to 
allow  the  leading  edge  to  defocus  slightly  for  the  sake  of  keeping  later 
points  in  better  focus. 
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5.  Error  of  Short  Impulse  Response  Approximation 

Landt,  Miller,  and  Van  Blaricum^  show  the  transient  echo  response  of  a 
thin  (hi-Q)  dipole.  Its  ringing  is  damped  after  5  cycles.  This  response  is  a 
good  example  to  gain  an  intuitive  idea  of  what  kind  of  range  is  needed  in 
the  approximation.  A  1/2  foot  dipole  should  ring  for  5  cycles  at  1  GHz. 
If  the  A/D  sampler  collects  data  at  2  Gs/s,  then  at  M-10,  all  of  the  5 
cycles  will  have  been  collected.  If  the  dipole  were  5  feet  long,  then  it 
would  ring  for  5  cycles  at  100  MHz.  So  M=100  to  collect  the  5  cycles. 
Generally,  high  frequencies  damp  quickly  and  low  frequencies  damp 
slowly. 

Figure  3  illustrate8  the  geometry  of  the  approximation.  The  discussion  will 
assume  this  geometry.  If  a  target  is  located  at  (i,k)=(237 ,0),  then  R  is  the 
distance  from  the  center  of  the  array  to  the  target,  al  is  the  distance  from 
the  end  of  the  array  to  the  target.  f237t(0)  is  the  perfectly-focused  data 
point  for  the  impulse  response  of  the  target.  The  next  point  (m=l)  in  the 
impulse  response  is  approximated  by  saying  that  al  ~  al  +  ad,  or  in 

general  a2~  al  +  maj.  If  e  is  the  difference  (error)  between  the 

approximation  and  the  actual,  then 
e  =  approximation  -  actual 

=  (al  +  mad)-«32  (9) 

-  Jr2  +:~-LsR  cos(0)  +  mad-J(R  +  mad)2+^--Ls(R  +  mad)cos(6) 


h.A.  Landt,  E.K.  Miller .  and  M.  Van  Blaricum  "WT -MBA/LLL1B:  A  COMPUTER  PROGRAM  FOR  THE  TIME- 
DOMAIN  ELECTROMAGNETIC  RESPONSE  OF  THIN -WIRE  STRUCTURES"  Lawrence  Livermore  Laboratory; 
May  6, 1974 
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Figure  3.  Geometry  for  calculation  of  approximation  error 

As  an  example,  suppose  R=2  Km,  6  -  90°,  Ls=l  Km,  and  m-10.  How 
many  degrees  off  (round  trip)  would  the  end  points  be  at  1  GHz? 
Plugging  into  equation  9,  we  find  £=.0447  meters.  So  the  round  trip 
phase  error  at  1  GHz  is  107°.  Figure  4  is  a  plot  of  the  error  as  a  function 

of  the  position  in  the  aperture  for  6=90,76,50,  and  30  degrees.  A  peculiar 
characteristic  is  that  the  error  peaks  at  76°.  This  peaking  is  a  result  of 
working  at  a  range  of  only  twice  the  aperture  length.  Figures  5  and  6 

show  the  error  as  a  function  of  the  beam  angle  6,  for  m  =  5,  10,  15,  and  20 
at  two  ranges. 
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Figure  6.  Error  as  a  function  of  6  with  R/Ls=6 


Slant  Range  =6000 
Sample  Rate=2GHz 

r=0 


Figure  7.  Error  at  right  aperture  end  as  a  function  of  m 


Slant  Range  =2000 
Sample  Rate=2GHz 
x=0 


Figure  7  is  a  plot  of  the  error  at  the  right  end  point  of  the  aperture  as  a 
function  of  m,  with  r=0.  The  error  at  the  right  end  is  greater  than  the 
error  at  the  left  end  for  0  angles  below  90°  -  so  this  is  the  worst  case. 
Figure  8  is  the  same  plot  but  with  v=na  and  n=4.  In  this  case,  the  leading 
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Degrees  off  a 1 1  GHz 


edge  of  the  impulse  response  (at  m= 0)  is  out  of  focus  by  about  45°  at  1 
GHz.  Perfect  focus  occurs  when  m=n=4. 


Figure  8.  Error  at  right  aperture  end  as  a  function  of  m 


Slant  Range  =2000 
Sample  Rate=2GHz 
T=4a 


To  conclude.  Figures  4  through  8  show  give  an  indication  of  the  bounds 
over  which  a  5 -cycle  ring  of  a  hi-  Q  scatterer  is  adequately  captured  by  the 
approximation  used  in  equation  4.  These  figures  also  show  that  making 
t  *  0  can  significantly  increase  the  depth  of  focus  on  resonant  targets. 
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6.  Efficient  Calculation 


Fast  focusing  can  be  broken  into  3  stages:  first,  pre-processing,  where 
interpolation  of  the  raw  data  is  performed;  second,  computing  a  set  of 
polynomial  coefficients  that  will  be  used  for  fast  index  calculation;  and 
third,  performing  the  focusing  summation  of  equation  7  using  the 
coefficients  found  in  second  stage  to  find  the  index  corresponding  to  the 
proper  time  shift. 

6.1.  Pre-Processing 

The  first  aspect  requiring  a  solution  is  a  method  to  efficiently  implement 
time  shifting.  The  SAR  data  is  collected  by  an  A/D  converter  that 
outputs  a  vector  of  N  numbers  (voltages)  at  each  aperture  position.  This 

vector  will  be  called  Sj(i ).  A  time  shift  is,  therefore,  simply  a  shift  in  the 

index  i  .  Typically,  the  time  shifting  obtained  by  indexing  on  the  original 
N-point  vector  is  not  fine  enough.  Finer  resolution  is  gained  by  a  two- 
stage  interpolator.  First,  a  high  quality  interpolator  is  used  to  produce  a 

new  K  point  vector  Sj(i).  It  is  usually  implemented  by  inserting  M-l 

zeros  between  each  data  point  in  the  original  vector  and  then  passing  the 
new  sequence  through  a  low-pass  FIR  (Finite  Impulse  Response)  filter. 
The  process  results  in  a  new  vector  with  nearly  the  desired  time- shift 
resolution.  The  length,  in  this  case,  is  K  =  MN .  Extra  fine  resolution  is 
gained  by  using  a  floating-point  index  with  simple  linear  interpolation  to 
find  a  value  between  any  two  data  points  in  the  interpolated  data  vector. 
For  example,  if  the  index  i  were  6.3,  then  the  interpolated  value  would 

just  be  .7sj(6)+.3Sj(7). 

The  focusing  algorithm  that  follows  uses  these  techniques  in  the  following 
sequence. 

1.  Read  in  one  N-point  vector  Sj(i )  of  raw  data. 

2.  Do  an  M-point  interpolation,  to  form  a  new  K-point  vector  Sj(i). 

3.  Iterate  for  all  pixels: 

a.  Find  the  floating  point  index  needed  for  a  pixel. 

b.  Find  the  data-value  for  that  pixel  by  either  rounding  the  index 
and  grabbing  the  value,  or  by  linearly  interpolating  between 
adjacent  values. 

c.  Sum  into  that  pixel  the  data-value  obtained. 
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6.2.  Fast  Index  Calculation 
6.2.1.  Approach 

Typically,  the  greatest  computational  load  in  backprojection  focusing  is  in 
finding  the  index.  Let  P(i,kj )  be  the  exact  index  needed.  Then  equation 
7,  the  2-D  focusing  equation,  becomes: 

fi,k  =  yLsj(P(i,krj))  (10) 

/ 

Calculating  P(i,k,j)  exactly  involves  3-D  trigonometric  solutions  for 
every  pixel  in  the  image  at  every  position  in  the  aperture.  Such 
calculation  is  prohibitive!  So  for  practical  implementation,  a  secondary 
approximation  is  utilized  to  speed  the  computations.  Extremely  efficient 
computational  methods  exist  for  finding  evenly  spaced  solutions  to 
polynomials.  Therefore,  the  approach  will  be  to  define  a  polynomial  in 

three  variables  G(i,kJ),  where  G(i,k,j )  =  P(i,k,j),  and  use  G  to  compute 
the  index. 

The  first  issue  that  arises  is  choosing  the  order  of  the  polynomial  needed 
for  each  variable.  For  now  (to  explain  the  method)  suppose  that  a  second 
degree  polynomial  is  adequate  for  each  of  the  3  variables  (i,k,j).  Now  the 
problem  can  be  re-stated  as,  "for  a  given  image  pixel  (i,k),  and  a  given 
aperture  position  (j),  find  the  coefficients  am  for  m=o..26  such  that 


G(i,kr  j)  ~  P(i,k,  j) 

=  %,k,j+%k,ji+cl2,k,ji2 

O  q  0  O 

=  [cotj+c\,jk+c2,jk  ]+[c3rj  +  c4/jk+c5'jk  ]i  +  [c6/j+c7/jk+c8/jk  ]t 

=  l(«0  +  #i  j+  #272)  +  (#3  +  a4)+  asjZ)k+  (a6  +  a7j  +a8j2)k2] 

+[(#9  +a10j+al  j/2 )  +  (fli2  +  «i3  7 + «14/2  >* + (*15  +  #167 + #17  i2  1* 

+[(#1  g  +  #19  j  +  #20/2  )  +  (#21  +  #22 7  +  #23/ 2  +  (#24  +  #25  /  +  #267  2  ^2 31'2 

Written  in  this  format,  it  is  easy  to  see  that  one  can  code  the  index 
calculation  as  loops  nested  3  deep;  with  j  as  the  outer  loop,  k  as  the  middle 
loop,  and  i  as  the  inner  loop. 

6.2.2.  Solving  For  The  Coefficients 

The  coefficients  am  can  be  found  numerically.  The  approach  is  to 
calculate  the  exact  index  required  for  M  points  in  the  image  at  each  of  H 
positions  in  the  aperture,  and  do  a  least-squares  fit  to  find  the  coefficients 
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G(i,k,j)  y= 


am  for  m=0..26.  A  generic  solution  to  this  problem,  for  arbitrary  image 
size,  is  found  by  using  a  pseudo  inverse  to  perform  the  least-squares  fit. 
When  matrix  A  is  not  square  but  rectangular,  then  A'1  is  known  as  the 
Pseudo  Inverse  and  is  defined  as: 

A-1  =(at  •  a)_1  •At.  (12) 

An  example  will  explain  the  method.  To  simplify  the  example,  we  will 
find  a  solution  for  only  a  single  position  in  the  aperture  at  j=> 0.  So  we  will 
find  the  cq  through  c8  needed  to  calculate  the  index  needed  for  so(G(i,k,0). 
Take  M=25  points  on  a  patch  to  be  focused:  5  ranges  (close-in  to  far-out: 
/=  0,  b,  2b,  3b,  4b)  on  each  of  5  azimuths  (left  side  to  right  side:  k= 0,  a, 
2a,  3 a,  4a).  The  vector  B  will  be  the  exact  (floating-point)  calculated 
index  needed  for  those  25  points  in  the  image.  The  matrix  A  will  be  set 
up  so  that  the  width  of  the  patch  is  4 a  and  the  depth  of  the  patch  is  4b.  So 
we  have 


Q~cO,Q+cl,ok  +  C2,okZ  +C3roi  +  c4:toik+c5rQik2  +  c6/0i2  +c7iQi2k  +  c8/0i2k2  (13) 


or 


G{i,k,j)  = 


1  k  k2  i  ik  ik 2  i2  i2k  i2k 2  *C, 


with 


A  «C  =  B 


(14) 


which  is. 
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Now  C  is  calculated  as  C  =  A_1B ,  and  is  used  to  find  a  floating-point 
index  pointing  into  the  data  array.  This  floating-point  index  is  either 
rounded  to  pick  the  nearest  point  or  is  used  to  do  a  linear  interpolation 
between  the  two  closest  data  points  in  order  to  find  a  data  value. 

Two  options  are  available  to  find  the  coefficients  oq— ^26-  The  first 
method  is  to  go  through  the  same  process  as  shown  but  for  the  full 

problem.  That  process  would  involve:  enlarging  C  to  hold  the  27  a0  26 
coefficients;  enlarging  B  to  hold  the  ideal  values  for  not  just  one  position 
in  the  aperture,  but  more,  say  5  positions;  and  enlarging  A  to  match. 

The  second  method  is  to  solve  for  vector  C  (as  shown)  at  a  number  of 
positions  in  the  aperture.  Then  construct  a  polynomial  in  j  for  each 
coefficient.  This  solution  does  not  involve  inverting  the  large  A  matrix 
but  will  result  in  a  less  than  optimum  solution. 

One  beauty  of  this  second  method,  is  that  it  lends  itself  to  correcting  for 
airplane  motion.  Since  A  is  independent  of  the  airplane  position,  it  is 
inverted  once.  The  matrix  multiplication  to  find  the  index  polynomial 
coefficients  can  be  applied  at  every  aperture  position  or  short  sub 
apertures  if  desired.  So  the  algorithm  allows  real-time  UWB  motion 
compensation. 

6.2.3.  Fast  Polynomial  Calculation 

Now  that  we  have  a  polynomial  to  find  the  index,  we  need  a  fast  way  to 
compute  the  polynomial.  Directly  calculating  an  N^1 2 3 4  degree  polynomial 
usually  requires  N  adds  and  N  multiplies.  If,  however,  a  sequence  of 
solutions  is  desired  with  the  variable  changed  in  fixed  increments  -  the 
case  we  have  here  -  then  more  efficient  means  are  available.  An 
algorithm  by  Nuttall9  describes  a  method  which  can  be  applied  here,  to 
calculate  the  index  recursively  without  the  multiplications.  The  procedure 
is  as  follows: 

1.  Let  X(i)  =  q0  +  qli+q2iz  be  the  polynomial  to  be  solved,  where 
;=0, 1,2,3.... 

2.  Let  JTj(0  =  X(0-X(/-1)  =  q,  ~q2  +2q2i 

3.  Observe  that  Aj  (/) — X1  (i  - 1)  =  2 q2 

4.  Therefore  a  recursion  can  be  set  up  where: 


5  Albert  H.  Nuttall,  "Efficient  Evaluation  of  Polynomials  and  Exponentials  of  Polynomials  for 
Equispaced  Arguments"  IEEE  ASSP-35  No.  10  (Oct.  1987)  ppl486-1487 
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X1  (i)  =  Xx  (i  - 1)  +  2 q2 ,  and  X(i)  =  Xx  (z  - 1)  +  X,  (i). 
The  starting  values  for  the  recursion  are: 


X(0)  -  qQ ,  and  Xl  (0)  -  qx  -  q2  +  lq2i  ^qx~q2 


The  same  derivation  can  be  applied  to  an  arbitrary  N1*1  degree  polynomial 
to  produce  a  recursive  formula  that  requires  N  adds  per  step.  The 
technique  can  also  be  expanded  to  the  multiple  dimensions  by  nesting. 
Appendix  B,  lists  all  routines  used  in  the  fast  focusing  algorithm. 
Subroutine  poly2  in  section  B.5  of  appendix  B  uses  this  technique 
directly. 

6.2.4.  Coefficient  Generator  Program 

Figure  9  is  a  diagram  showing  how  the  image  is  broken  down  along  with 
the  names  used  in  the  subroutines. 


k  boxes=3 


i  boxes=2 


oe 

6  oxb  d>  oA  d>  a4> 

hOX)OOX)QO% 

e-©G-e-©GH& 

T  0.00  0.0  0  0, 
ooboooq) 
oeo&e-o 


k_pixes=3 


_pixes=4 


Figure  9.  Image  Partitioning  and  Notation 


The  basic  flow  chart  is  shown  in  figure  10.  Each  time  the  inner  loop  goes 
through  all  the  n’s.  The  Zp's  are  the  aO  to  a26  coefficients  for  the  box 
(defined  by  h  and  m)  and  the  sub-aperture  (defined  by  1).  The  flow-chart 
shows  that  the  polynomial  degree  for  the  Cn  is  not  always  the  same.  The 
degree  was  optimized  based  on  the  errors  caused  by  going  to  a  lower 
degree.  Appendix  A  provides  a  complete  listing  of  the  program  that 
generates  the  focusing  coefficients  from  an  input  file  with  the  geometry. 
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Figure  10.  Coefficient  Generator  Flow  Chart 

6.3.  Fast  Focusing  Algorithm 

The  fast  focusing  algorithm  simply  applies  the  fast  polynomial  solution 
technique  to  find  the  indexing  and  does  the  signal  summation  into  all  the 
pixels  in  the  image.  As  outlined  above,  the  index  calculation  is  done  in 
three  nested  loops.  To  aid  in  explanation,  subroutines  Inner_Loop, 

Middle . Loop,  and  Outer_Loop  are  pseudo-coded10  examples  of  how  the 

focusing  routine  is  written  using  the  recursive  technique  with  nesting.  A 
second  degree  polynomial  is  used  throughout  for  illustration  purposes. 
These  subroutines  assume  that  the  image  is  broken  into  a  number  of  boxes 
where  each  box  has  its  own  set  of  coefficients.  lnner_Loop  is  the  simplest 
subroutine  and  follows  the  derivation  of  the  recursive  formula  directly.  It 
essentially  increments  i  to  focus  a  single  line  in  a  box. 


Inner  Loop  Subroutine 

float  *f; 

pointer  to  first  pixel  in  a  line;  f(Lstart,k) 

float  *s; 

pointer  to  signal-data  vector  from  jth  aperture 
position 
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float  * 


ointer  to  coefficients  vector  for  index  calculation 


int  s_max; 

s  range  is  s!01..sls_maxl  s„max; 

int  Lpixes; 

size  of  the  box  (pixels)  along  i  axis  is  Lpixes; 

f  stop  i  ~f  +  i  pixes-1; 

T=q[2]+q[2]; 

Rl=q[l]-q[2]; 

R0=q[0J; 

Set  up  initial  conditions 

If  R0  <  -0.5  then  repeat 
{f++; 

R1=R1+T; 

R0=R0+R1}; 
until  R0> -0.5}; 

If  (f  >  f_stop  i)  then  return; 

Perform  Clipping  when  indexcO 
(prior  to  beginning  of  actual  data) 

repeat 


f++; 


R1=R1+T; 


index=Round(RO); 
if  index>s.max  return; 


*f  =  *f  +;*(s+index' 


until  (f  =  fjstop); 
return; 


Perform  1st  summation 


begin  loop  for  a  line  in  box 


increment  pointer  to  next  pixel 


increment  index  R0=R0+R1; 


(clip  when  index  is  beyond  actual  data  record 
length) 


sum  into  pixel  the  new  data; 


(END  InnerJLoop) 


Nesting  of  the  recursive  approach  means  that  the  coefficients  qO,  ql,  and 
q2  are  each  formed  by  polynomials  in  k.  Each  polynomial  is  incremented 
recursively  in  the  subroutine  Middle_Loop.  So  Middle_Loop  sums  into  a 
single  box  the  data  from  a  single  aperture  position. 


Middle_Loop(f,  s,  c,  s_max, 
Lpixes,  kjnc,  box . offset); 

Subroutine 

float  *f; 

pointer  to  first  pixel  of  box 

float  *s; 

pointer  to  jth  data  vector 

float  *c; 

pointer  to  coefficients  vector  for  index  calculation 

int  s„max; 

s  range  is  s[0]..s[s_max] 

int  k_inc; 

F[i,k]  =  *(f  +  i  +  k  *  k_inc)  so  k . inc  is  the  total  range  line 

length 

int  Lpixes; 

size  of  the  box  (pixels)  along  i  axis 

int  box_offset; 

Number  to  add  to  f,  to  move  pointer  to  last  line  in  box 

Llast  =  f  +  box_offset; 

f_last  =  address  of  the  first  pixel  in  the  last  line  of  box 

qO  T=c[2]+c[2]; 
qO  Rl=c[l]-c[2]; 
q[Q]rcLQI; _ 

Set  up  initial  conditions  for  qO 
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Set  up  initial  conditions  for  ql 

q2  T=c[8]+c[8J; 

q2_Rl=c[7]-c[8]; 

qr21=cr6i; 

Set  up  initial  conditions  for  q2 

Call  Inner JLoop(f,s,q..); 

Perform  summation  on  first  line 

repeat 

start  loop  for  rest  of  lines 

f  =  f  +  k_inc; 

increment  pointer  to  next  line  in  box 

qO  Rl=qO  Rl+qO  T; 
q[0]=q[0]+q0_Rl; 

iterate  qO 

ql  Rl=ql  Rl+ql  T; 
q[l]=q[l]+ql_Rl 

iterate  ql 

q2  Rl=q2  Rl+q2  T; 
q[2]=q[2]+q2_Rl; 

iterate  q2 

Call  Inner_Loop(..); 

Perform  summation  on  current  line  in  box 

Until  (f  =  f Jast); 

Return; 

Loop  until  all  lines  done 

Again,  nesting  of  the  recursive  approach  means  that  the  coefficients  c0..c8 
are  each  formed  by  polynomials  in  j.  Each  polynomial  is  incremented 
recursively  in  the  subroutine  OuterJLoop.  So  OuterJLoop  sums  into  a 
each  box  the  data  from  a  multiple  aperture  positions.  Since  it  is  desirable 
to  only  read  the  signal-data  once,  OuterJLoop  will  both  call  Middle_Loop 
for  each  box  and  increment  the  coefficients  for  each  box  separately  as  it 
increments  j  across  the  aperture. 


Outer  Loop(f,a,s_max,k_boxes,i_boxes,k_pixes,i_pixes,i_beg,  Lend) 

float  *fo; 

pointer  to  first  pixel  in  entire  image;  all 
pixels  set  to  zero. 

float  *a; 

pointer  to  array  of  coefficient  vectors  for 
index  calculation;  one  vector  per  box. 

int  sjmax; 

s  has  range  of  sfO]  to  sf s_max] 

int 

number  of  boxes  in  k  axis  k_boxes; 

int 

number  of  boxes  in  i  axis  i_boxes; 

intk  pixes; 

size  of  the  box  (pixels)  in  k  axis 

int 

size  of  the  box  (pixels)  in  i  axis  i_pixes; 

int  i  beg; 

start  of  data  vectors  to  be  processed 

int)  end; 

stop  of  data  vectors  to  be  processed 

k Jnc  =  i_pixes  *  i_boxes; 

(F(i,k))=  *(f+i+k*kjnc) 

box  offset  =  k  inc  *  (k_pixes-l); 

Used  by  MiddleJLoop 

box  inc  k  =  Lpixes  +  box_offset; 

box  Jnc  J  =  Lpixes; 

j~j_start; 

read  sf0..s_.maxj  for  i  th  position; 

read  signal  data  for  first  aperture  position 
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f  =  fo  -  (box_Jnc  Jc+box_inc_i); 
box=-2; 

for  k_box  =  0  to  k_boxes-l; 
box  =  box  +  1 ; 
f  =  f  +  box_inc_k; 
for  i_box  =  0  to  i_boxes-l; 
box  =  box  +  1; 
f  =  f  +  box  inc  i; 


cO_T[box]=a[2,box] +a(2,box] ; 
cO_R  1  [box]=af  1  ,box]-a[2,box] ; 
c[0,box]=a[0,boxl; 


c  1_T[  box] =a[5 ,  box]  ^a[5 ,  box] ; 
clJRl[box]=a[4,box]-a[5,box]; 
cl  c[l,box]=a[3,box]; 


set  up  to  loop  through  all  boxes 


c8_T[box]=a[27,box]+a[27,box]; 
c8_R  1  [box]=a[26,box]-a[27,box]; 
c  [8 ,  box] =a[25 ,  box] ; 


Call  Middle  Loop(f,s,c[box 


repeat 


+1; 


read  s[0..s_max]; 


f  =  fo  -  (box  inc  k+box  inc  i 


box=-2; 

for  k_box=0  to  k_boxes-l; 
box  =  box  +  1; 
f  =  f  +  box_inc_k; 
for  i_box=0  to  i„boxes-l; 
box  =  box  +  1 ; 
f  =  f  +  box  inc  i; 


cO_Rl  [box]=cO_Rl  [box]+cO_T[box]; 
c[0,box]=c[0,box]+c0__R  1  [box]; 


cl_Rl[box]=cl_Rl[box]+cl_T[box]; 
c  [  1  ,box]=c[  1  ,box]  +cl_Rl  [box] ; 


c8_R  1  [box]=c8_R  1  [box]+c8_T[box]; 
c  [8  ,box] =c[8,  box]+c8_R  1  [box]; 


initial  cond.  for  cO  for  current  box 


initial  cond.  for  current  box 


initial  cond.  for  ,c2..c7  for  current 
box 


initial  cond.  for  c8  for  current  box 


first  summation  for  current  box 


loop  through  all  boxes 


loop  through  all  aperture  positions 


increment  aperture 


position  read  data  for  that  aperture 
osition 


set  up  to  loop  through  all  boxes 


iterate  cO  for  current  box 


iterate  cl  for  current  box 


iterate  c2..c7  for  current  box 


iterate  c8  for  current  box 


Sum  into  current  box 


Loop  until  all  boxes  done 
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Until  j=j__stop; 

Loop  until  all  aperture  positions  are  done. 

Return  (END  Middle_Loop) 

6.3.1.  Polynomial  Order  Selection  And  Computational  Load 

A  study  was  conducted  to  determine  the  effect  of  polynomial  order  on  the 
image  size  and  aperture  length.  Referring  to  Figure  3,  the  scenario 

geometry  was  for  an  aperture  of  Ls= 384  feet,  a  squint  angle  of  0  =  85°,  a 
slant  range  of  550  feet  to  the  closest  range  line  (=0,  an  elevation  of  60 
feet,  an  image  area  54°  wide,  for  the  worst-case  position  and  box  j~ 3  and 
(i,k)=( 0,3),  and  for  a  4  Gs/s  data  rate 


with  record  lengths  (s . max)  of  4096  points.  If  the  error  in  the  floating 

point  index  is  restricted  to  ±.  5 ,  then  this  geometry  produces  the  following 
results: 


Poly  Nom 
Degree 

max 

i_pixes 

max 

k_pixes 

1 

228 

27 

2 

964 

105 

3 

2191 

315 

4 

657 

By  using  the  nested  recursive  technique  as  shown  in  the  example 
subroutines,  the  computational  load  is  computed  as  shown  below. 


i’ 

polynomial  order  for  i  index; 

k’ 

polynomial  order  for  k  index; 

1' 

polynomial  order  for  i  index; 

J 

Number  of  points  in  the  aperture; 

I 

Number  of  range  bins  (pixels)  in  box  (Lpixes); 

K 

Number  of  azimuth  bins  (pixels)  in  box  (k_pixes); 

B 

Number  of  boxes 

L0=I*(i’+l) 

Adds  per  Inner  Loop  call 

Ll=K*(Lo+(i'+l)k’) 

Adds  per  Middle  Loop  call 

L=J*(Ll+(i’+l)(k’+l)f) 

Adds  per  Outer  Loop  (adds  per  image) 

Grouping  terms  we  get: 

L=(inner  loop  adds)  +  (middle  loop  adds)  +  (outer  loop  adds) 

=IKJB(i'+l)  +  KJB(if+l)k'  +  JB(i'+l)(k'+l)j' 

From  this  equation  it  is  clear  that  a  low  order  polynomial  is  crucial  in  the 
inner  loop.  Conversely,  the  order  of  the  polynomial  in  the  outer  loop  can 
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be  high  with  little  impact  on  the  overall  speed.  The  table  below  shows  the 
computational  load  for  various  configurations  with  j'=3  and  J=2304.  The 
upper  number  is  the  tonal  number  of  adds  (L)in  Giga-adds.  The  lower 
number  is  the  total  number  of  boxes  in  the  image  (B).  The  center  number 
is  the  bytes  of  storage  needed  for  the  coefficient  tables  which  is 
(i'+l)(k’+l)(j'+l)*B*(4  bytes/word). 


Table  2.  Computational  Load  Versus  Partitioning 


i’ 

1 

2 

3 

4 

K 

27 

103 

256 

512 

boxes 

across 

19 

5 

2 

1 

total  pixels 
across 

513 

515 

512 

512 

q 

I 

boxes 

down 

total  pixels 
down 

i 

227 

18 

4086 

L=9.7 1 1 

S=21888 

B=342 

L=9.786 

$=8640 

B=90 

L=9.769 

S=4608 

B=36 

L=9.810 

S=2880 

B=18 

2 

817 

5 

4085 

L=14.51 

S=9120 

B=95 

L=14.58 

$=3600 

B=25 

L=14.51 

S=1920 

B=10 

L=14.53 

S=1200 

B=5 

3 

2043 

2 

4086 

L=19.33 

S=4864 

B=38 

L=19.41 

S=1920 

B=10 

L=19.31 

S=1024 

B=4 

L=19.32 

S=640 

B=2 

6.4.  Implementation  Notes 

The  algorithm  developed  here  requires  no  multiplication  except  in  the  pre¬ 
processing  (interpolation)  stage.  There  already  exist  DSP  chips  whose 
architecture  is  ideal  for  the  FIR  interpolation  task.  The  bulk  of  the 
computations,  however,  occur  in  the  summing  and  index  calculation  stage. 
Two  facts  make  the  summing  algorithm  attractive.  First,  adders  take 
considerably  less  area  to  implement  in  a  VLSI  chip  than  multipliers.  And 
second,  parallel  operation  is  extremely  simple;  the  image  can  simply  be 
broken  into  boxes  with  a  separate  processor  working  independently  on 
each  box.  Except  for  passing  the  signal-data  to  each  processor,  the 
processors  could  run  independently;  only  the  coefficient  table  would  be 
different.  Therefore,  design  of  a  custom  LSI  chip  appears  practical  and 
could  result  in  real-time  speeds. 

If  an  off-the-shelf  DSP  chip  with  a  parallel  adder  and  multiplier  is  used  to 
perform  this  algorithm,  then  other  options  become  available.  For 
example,  the  multiplier  can  be  used  in  the  summing  algorithm  to  do 
interpolation  on  the  indexing  rather  than  rounding.  The  multiplication 
required  can  be  done  during  other  add  cycles  such  that  it  takes  zero  time. 
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Another  possibility  is  to  calculate  the  index  polynomial  directly  instead  of 
recursively.  Finally,  all  the  operations  can  be  fixed  point.  So  the  address- 
generator  ALU  (Arithmetic  Logic  Unit)  can  sometimes  be  used  in  parallel 
with  the  floating  point  units. 
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7.  Beam  Patterns  and  Side  Lobe  Structure 


To  develop  an  intuitive  understanding  for  the  beam  shape,  16  figures  were 
assembled  to  display,  in  the  time  domain,  the  main-lobe  and  side-lobe 
beam  patterns.  Simulation  of  a  resonant  target  was  done  by  generating  a 
data  record  for  every  position  in  the  aperture.  The  simulated  scene  was  a 
single  point  target.  The  target  bearing,  refering  to  figure  3,  was  squinted 
15°  off  broadside  (0=105°).  Range  to  the  target  R  was  750  feet.  The 
aperture  length  Ls  was  385  feet.  The  aperture  height  was  60  feet.  The 
point  target  had  an  impulse  response  of 


s(t)  ■ 


sin(2^ft)[.5+.5cos(4^it)]  for0<r< 

<-U)°.23/  , 

sin(2n/i)e  forr>  — 


4/ 


(15) 


The  data  simulated  a  2  GHz  sample  rate.  A  record  length  of  2048  samples 
was  made  for  each  aperture  position.  Each  record  was  pre-processed  with 
a  times- 8  interpolator  to  produce  16K  records.  Interpolation  was  done  by 
the  standard  FIR  (Finite  Impulse  Response)  filter  method.  A  255  tap 
Parks-McClellan  low-pass  filter  with  a  950  MHz  cutoff  was  used  to  do  the 
interpolation.  Equation  7  (with  n= 0)  was  used  to  produce  the  focused 
beams.  A  Hilbert  Transform  was  used  to  obtain  the  magnitude  that  is 
plotted  in  the  figures.  Plots  were  made  on  a  dB  scale.  Sixteen  plots  were 
made  to  fill  the  matrix  shown  in  table  3. 


Table  3.  Point  Response  Function  Plots 


Frequency 

On-Axis 

On-Axis 

Above-Axis 

Above-Axis 

Equal 

Equal 

Hamming 

Equal 

Hamming 

50MHz 

Figure  11 

Figure  12 

Figure  13 

Figure  14 

200MHz 

Figure  15 

Figure  16 

Figure  17 

Figure  18 

400MHz 

Figure  19 

Figure  20 

Figure  21 

Figure  22 

900MHz 

Figure  23 

Figure  24 

Figure  25 

Figure  26 

The  on  axis  plots  were  made  so  that  amplitude  values  could  be  easily  read 
off  of  the  plots.  The  above  axis  plots  were  made  in  order  to  see  the 
sidelobe  structure,  and  the  ring-down  time  of  the  target.  The  axis  labels 
were  shifted  so  that  0  range  was  at  the  point  target  and  so  that  0  degrees 
was  the  bearing  centered  on  the  target.  Notice  that  all  plots  are  3-D,  even 
the  on-axis  plots.  The  various  horizontal  lines  in  the  on-axis  plots  are  the 
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beam  pattern  on  the  range  bin.  These  horizontal  lines  are  identical  to 
those  shown  in  the  above  axis  plots;  they  are  just  being  viewed  "end  on." 

These  figures  are  unique  in  that  they  show  how  the  sidelobes  spread  in 
time  as  one  moves  off  of  the  main  beam.  Due  to  the  shaip  rise-time  of  the 
target  echo,  the  contribution  from  the  near  and  far  ends  of  the  aperture  can 
be  seen  as  the  early  and  late  peaks  in  the  sidelobe  structure.  Hamming 
weighting  was  applied  to  the  array  to  enable  one  to  see  the  effect  of 
weighting  on  the  sidelobe  structure.  It  is  interesting  to  note  that  although 
the  peak  sidelobe  levels  only  drop  marginally,  the  average,  or  integrated 
sidelobe  levels  are  significantly  lower  when  tapered  aperture  weighting  is 
used.  The  weighting  also  reduces  the  near  and  far  peaks  in  the  sidelobes 
since  the  weighting  reduces  the  contribution  of  the  array  ends. 

Since  the  antenna  is  a  linear  system,  it  is  no  surprise  that  the  beam  pattern 
behaves  basically  as  classic  antenna  theory  predicts.  The  aperture  beam 
pattern  is  a  function  of  1)  how  long  the  aperture  is  relative  to  wavelength, 
and  2)  what  kind  of  weighting  is  used  to  sum  the  points  in  the  aperture.  If 
uniform  aperture  weighting  is  used  (a  straight  summation),  then  the 
typical  sin(/tyr)/(/fyO  pattern  occurs,  where:  \j/  is  the  angle  off  the  main 
beam;  and  «=  Ls  /  X .  establishes  the  width  of  the  main  beam.  As  the 
wavelength  gets  small,  or  as  the  aperture  gets  long,  the  beam  gets  narrow. 
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Image  Amplitude  (dB)  Image  Amplitude  (dB) 


2  1.5  1.0  .5  0  .5  1.0  -1.5  -2.0 

Azimuth  (degree) 


Figure  11. 

50MHz  Point  reflector 
with  Equal  Weighting 
End  View 


'2  1.5  1.0  .5  0  .5  1.0  -1.5  -2.0 

Azimuth  (degree) 


Figure  12. 

50MHz  Point  reflector, 
Hamming  Weighting, 
End  View 
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50MHz  Point  reflector. 
Equal  Weighting, 
Above  Axis 


Figure  14. 

50MHz  Point  reflector 
Hamming  Weighting 
Above  Axis 


Image  Amplitude  (dB)  Image  Amplitude  (dB) 


Figure  15. 

200MHz  Point  reflector 
Equal  Weighting 
End  View 
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Figure  16. 

200MHz  Point  reflector 
Hamming  Weighting 
End  View 
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Image  Amplitude  (dB) 


© 


Image  Amplitude  (dB)  Image  Amplitude  (dB) 
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Azimuth  (degree) 


Figure  19. 

400MHz  Point  reflector 
Hamming  Weighting 
End  View 


© 


o 

cst 


o 


© 

■*r 


'Si'i 

1 


4=  Ml 

1*1  'J  fijj 


Ilifi,:," 


Jf 


'fi’rti;; 

l 


'm  . 


1 . 5 


1 . 0 


.5  0  .5  1.0 

Azimuth  (degree) 


1.5  -2.0 


Figure  20.  ' 

400MHz  Point  reflector 
Hamming  Weighting 
End  View 
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Figure  21. 

400MHz  Point  reflector 
Equal  Weighting 
Above  Axis 


Figure  22. 

400MHz  Point  reflector 
Hamming  Weighting 
Above  Axis 


Figure  23. 

900MHz  Point  reflector 
Equal  Weighting 
End  View 
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Figure  24. 

900MHz  Point  reflector 
Hamming  Weighting 
End  View 
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Image  Amplitude  (dB)  -  Image  Amplitude  (dB) 


Figure  25. 

900MHz  Point  reflector 
Equal  Weighting 
Above  Axis 


O 


O 


Figure  26. 

900MHz  Point  reflector 
Hamming  Weighting 
Above  Axis 
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8.  Conclusion 


This  report  has  identified  a  frequency  independent  processing  algorithm  to 
get  a  perfectly  focused  impulse  response  from  any  object  at  any  position 
from  an  ultra  wide  bandwidth  synthetic  aperture  radar.  The  report  also 
presented  an  approximation  that  reduces  the  computational  complexity  of 
the  algorithm.  An  error  analysis  of  this  approximation  demonstrated  its 
applicability  to  focus  even  hi-Q  objects  in  the  near  field  of  an  aperture. 
An  implementation  that  is  both  computationally  efficient  and  applicable  to 
real-time  motion  compensation  was  also  presented.  Finally,  plots  were 
made  demonstrating  the  capability  of  the  algorithm  to  focus  ringing 
targets  over  an  18-to-l  bandwidth. 
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APPENDIX  A 


Coefficient  Generator  Program 

A.l.  Main  Program  to  Calculate  Coefficients  (main.c) 


Program :  coefgen 

Description  :  This  is  the  main  program  to  generate  coefficients 
for  fast  focusing  algorithm  . 


ifc*************^*************************************************************/ 


#include  <malloc.h> 

#include  <stdio.h> 
include  <math.h> 

#include  "coefgen. const" 

#include  "coefgen.var" 

main() 

{ 

long  cacheset; 

FILE  *fp,*tstart  fp; 
char  fn[80]; 
long  i; 

ia=(BOX_SIZE„AZIMUTH-4)/4; 

ib=(BOX_SIZE_RADIAL-3)/4; 

/*  initialize  variables  */ 
coefgen_varinit(); 

/*  initialize  pseudo  inver  matrix  */ 
initmat(); 


printf("Enter  coefficient  a  output  file  name  ===>  "); 

scanf("%s",fn); 

fp=fopen(fn,"w"); 

if(fp==NULL) 

{ 

printfC’Error  open  output  file  \n"); 
exit(l); 

} 

/* 

printf("Enter  data  acquisition  timing  file  name  =>  "); 
scanf("%s",fn); 

*/ 
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/* 

tstart_fp=fopen("delaytime.dat’',"r"); 

if(tstart_fp=-NULL) 

{ 

printf( "Error  open  input  file  \n"); 
exit(O); 

} 

read_delaytime(tstart_fp); 

*/ 

coefgen(fp); 

} 


A.2.  Main  Coefficient  Generator  Routine  (coefgen.c) 

#include  "coefgen.h"  /*geometry  file*/ 

#include  "coefgen.var"  /*list  of  variables*/ 

#include  <stdio.h> 

#include  <math.h> 
extern  indexQ; 
extern  mxmul_(); 
extern  poly_fit(); 

/*  Note  that  Reference  point  is  now  taken  at  coordinate  kpixel=  n_azimuth/2- 1  */ 
/*  ipixel=  n_radial/2  -1  */ 

/*  At  every  position  data  were  taken  so  that  the  reference  point  is  at  center  */ 

/*  data  buffer  */ 

void  coefgen(fp) 

FILE  *fp; 

{ 

float  input[NPOINT_APER/4]; 
long  i,j; 

long  section, kbox,ibox,position,position_start,position_stop; 
long  ipoint,kpoint; 

for(i=0;i<NPOINT„APER/n„section;i++) 

{ 

input[i]=(float)i; 

} 

for(section=0;section<n_section;section++) 

{ 

position_start=section*n_aper/n_section; 

position . stop=position . start+n_aper/n_section; 

for(kbox=0;kbox<NBOX_AZIMUTH;kbox++) 

{ 

for(ibox=0;ibox<NBOX_RADIAL;ibox++) 

{ 

printfC’processing  section  %d  ibox  %d  kbox  %d\n",section,ibox,kbox); 
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for(position=position_start;position<position„stop;position++) 

{ 

/******CalcuIating  Actual  index  for  all  samples  points  in  a  box**********/ 
i=0; 

for(ipoint=0;ipoint<NPOINT_BOXSAMPLE_RADlAL;ipoint++) 

{ 

for(kpoint=0;kpoint<NPOINT_BOXSAMPLE_AZIMUTH;kpoint++) 

{ 

ipixel=ibox*BOX_SIZE_RADIAL+ipoint*ib; 

jpixel=kbox*BOX_SIZE_AZIMUTH+kpoint*ia; 

d=((float)position-((float)n_aper/2.-l.))*(aper_length/((float)n__aper-l.));  /*  distance  from 
radar  to  center  of  aperture  */ 
x=sqrt(rcenter_ref*rcenter_ref-radar_height2); 

r_ref=sqrt(radar_height2+x*x*c_theta_center_ref*c_theta__center_ref-f(d- 

x*s . theta_center_ref)*(d-x*sjheta_center_ref)); 

rmin=r_ref-d_rcenter*((float)n_radial/2.- 1 .); 

index(&ipixel,&jpixel,&findex[i]); 

i++; 

} 

} 

/*  generate  coef  for  this  box  at  this  position  */ 

mxmuls(mat,&mat_c_stride,&stride,findex,&findex_c_stride,&stride,coef,&coef_c_stride,&stri 

de,&n_coef,&i_one,&n_boxsample); 

/*  save  in  coefc[i]  [position]  */ 

for(i=0;i<COEF_SIZE;i++) 

{ 

coefc[i]  [position-position_start]  =coef[i] ; 

} 

/*  for  debug  only  */ 

/*  printf("pos-  %d  tstart=  %.3e\n", position, 2.0*rmin/c);  */ 

/*  end  debug  */ 

}  /*  position  */ 
for(i=0;i<COEF_SIZE;i++) 

{ 

printf("Doing  polyfit  for  coef  %d\n",i); 
poly_fit(input,&coefc[i3[0],n_aper/n_section,deg[i],&coefa[i][0]); 
for(j=0;j<(deg[i]+l);j++) 

{ 

fprintf(fp,''%e\n",coefa[i][]]);  /*  save  a[i]  to  file  */ 

} 

} 
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}  /*  ibox  */ 
}  /*  kbox  */ 
}  /*  section  */ 


}  /*  subroutine  */ 


A.3.  Subroutine  Find  Ideal  Index  Vector  (index.c) 


Subroutine: 

Input: 


index.c 


ipixel 

pixel  radial  coordinate 

jpixel 

pixel  angle  coordinate 

a__span 

angle  spanned  by  the  patch 

a_ofset 

offset  angle  of  the  center  line 

d 

distance  from  middle  position 

d_rcenter 

sampling  range 

dr8 

(sampling  range)/8 

rmincenter 

min  range  rom  center  position 

hgt 

height  of  building 

hgt2 

square  the  height 

length 

>ut: 

length  of  the  building 

findex 

floating  point  index  to  ivalue  of  that  pixel 

#include  "coefgen.const" 
#include  "coefgen.var" 
#include  <math.h> 

index(ipixel,kpixel,fmdex) 
int  *ipixel,*kpixel; 
float  *findex; 

{ 


theta=a_of set-a_span/2. +(float)  *kpixel  *d_theta; 

c_theta=cos(theta) ; 

s_theta=sin(theta); 

rcenter=rmincenter+(float)  *  ipixel  *d_rcenter; 
x=sqrt(rcenter*rcenter-radar_height2); 

r=sqrt(radar„height2+x*x  *c_theta*c_the  ta+(d-x*s_theta)  *  (d-x  *  s_theta)); 
sicfmdex=(r-nnin)/drB; 

} 
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A.4.  Geometry  File  -  Declare  All  Global  Constants  (coefgen.h) 


^define  N POI NT_  AZIMUTH  600  /*  number  of  bearing  lines  */ 

#define  NPOINT_RADIAL  4095  /*  number  of  radial  lines  */ 

#defme  NPOINT_APER  2304  /*  #  of  positions  in  aperture  */ 

#define  NPOINT_DATA  2048  /*  number  of  original  data  points  */ 

#deflne  NPOINT_D A TA_INTER  NPOINT_DATA*8  /*  number  of  data  points  after 

interpolated  */ 

#defme  PI  3.141592654 

#define  RADAR_HEIGHT  60.  /*  radar  height  in  feet  */ 

#define  A_OFSET  -15.0  /*  offset  angle  from  center  line  */ 

#define  A_SPAN  54.0  /*  spanning  angle  of  the  patch  */ 

#define  RMINCENTER  550.  /*  range  from  center  position  to 

nearest  point  on  the  patch  */ 

#define  RCENTER_REF  802.0  /*  range  (ft)  of  ref.  point  from  center  pos.  */ 

#define  THETA_CENTER„REF  -15.0  /*  angle  (deg)of  ref.  point  from  center  pos.  */ 
#define  SAMPLING  RATE  2.0e09  /*  sampling  frequency  of  signal  */ 

#defme  SAMPLING.PERIOD  5.0e-10  /*  ts=l/fs  */ 


#define  BOX„„SIZE_RADIAL  195 
#define  BOX_SIZE_ AZIMUTH  100 

#define  NBOX  RADIAL  NPOINT_RADIAL/BOX_SIZE_RADIAL 
#define  NBOX_AZIMUTH  NPOINl'_AZlMUTH/BOX_SIZE_AZLMUTH 

#define  SECTION_SIZE  4  /*  #  of  sections  for  one  aperture  */ 

#define  COEF_SIZE  6  /*  #  of  coeffs  for  curve  fit  */ 

#define  MAXDEG  3  /*  maximum  degree  for  poly,  fit  */ 

Mefine  NPOINT_BOXSAMPLE_RADIAL  5  /*  #  of  samples  in  a  box  in  radial 

direction  */ 

#define  NPOINT_BOXSAMPLE_ AZIMUTH  5  /*  #  of  samples  in  azimuth  direction 

for  every  box  */ 

#define  NPOINT_BOXSAMPLE 

NPOINT_BOXSAMPLE_RADIAL*NPOINT„BOXSAMPLE_AZIMUTH 

/*  #  of  samples  in  a  box  for  curve  fit  */ 

A.5.  Declare  All  Global  Variables  (coefgen.var) 

struct  complex  {  float  r,i; } ;  I*  define  a  complex  type  *1 

/*  Data  Variables  */ 

float  pixel[NPOINTAZIMUTH]  [NPOINTJRADIAL] ;  /*  array  of  image  */ 

float  data[NPOINT_DATA];  /*  original  data  bufer  */ 

float  data_inter[NPOINT_DATA„INTER];  f*  interpolated  data  buffer  */ 

struct  complex  data Jnter_fft[N  POINTED  AT  A_INTER/2];  /*  Real->Complex  Forward  FFT 
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of  datajnter  *i 

float  filter_coef[NPOINT_p  AT A_INTER] ;  /*  filter  coefficients  */ 

struct  complex  filter_fft[NPOINT  DATA_INTER/2j;  /*  Real->Complex  Forward  FFT 

of  filter_coef  */ 

float  nyquist_filter;  /*  value  of  FFT  of  filter  at  Nyquist  point  */ 
float  nyquist_data;  /*  value  of  FFT  of  data  at  Nyquist  Point  */ 


/*  Geometry  Variables  */ 

float  a_ofset;  /*  offset  angle  from  the  center  line  */ 

float  a_span;  /*  spanning  angle  of  the  patch  */ 

float  d_theta;  /*  delta  angle  */ 

float  x, theta;  /*  coordinate  of  a  pixel  */ 

float  c_theta,s_theta;  /*  cosine  and  sine  of  theha  */ 


float  d;  /*  distance  from  radar  to  center  position 

positive  to  the  right,  neg.  to  the  left  */ 
float  reenter ,rmincenter,rmaxcenter;  /*  range  from  center  position  */ 
float  rcenterjref;  /*  range  from  center  pos.  to  ref  point  */ 

float  theta_center_ref;  /*  angle  of  ref.  point  from  center  position  */ 


float  r_ref;  /*  range  from  any  pos.  to  ref.  point  */ 

float  d_rcenter;  /*  sampling  distance  from  center  position  */ 

float  dr8;  (*  (sampling  distance)/8  */ 

float  r,rmin,rmax;  /*  range  from  any  position  */ 

float  aperjength;  /*  length  of  aperture  */ 

float  radar_height;  /*  height  of  radar  */ 

float  radar_height2;  /*  square  the  height  of  radar  */ 

/*  Radar  Variables  *1 

float  c;  /*  speed  of  wave  */ 

float  ts;  /*  sampling  period  of  signal  */ 

float  fs;  /*  sampling  frequency  of  signal  */ 


long  n_azimuth; 
long  njradial; 
long  n_aper; 
long  pix_size; 
long  n_data; 
long  n_data_inter; 
long  n_data_inter_half; 
long  n_boxsample; 
long  n_coef; 
long  n_section; 


/*  #  of  points  in  azimuth  direction  */ 
/*  #  of  points  in  radial  direction  */ 

/*  #  of  positions  in  aperture  */ 

/*  #  of  points  in  pixel  array  */ 

/*  number  of  original  data  points  */ 

/*  number  of  interpolated  data  points  */ 
/*  1/2  #  of  interpolated  data  points  */ 
/*  #  of  samples  in  a  box  for  curve  fit  */ 
/*  #  of  coefficients  curve  fit  */ 

/*  #  of  sections  for  one  aperture  */ 


/*  SSL  VAriables  */ 
float  f_zero; 
float  f_one; 
long  i_one; 


/*  floating  point  zero 
/*  floating  point  1 
/*  integer  1 


*/ 

*/ 

*/ 
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long  stride; 


/*  stride  for  supercard  ssl 


*/ 


/*  working  variables  */ 
long  n_position; 

float  e_index,a_index,  error,  maxerror.minerror; 
long  ierr_max,jerr_max,ierr_min,jerr__min; 
long  column_max,row_max,column_min,row_min; 
long  ia,ib; 
long  ipixel  jpixel; 

float  c_theta_centerjre  f,  s_the  ta_cen  ter . ref ; 

float 

mat[COEF_SIZE][NPOINT_BOXSAMPLE],findex[NPOINT_BOXSAMPLE],coef[COEF_SIZ 

E]; 

float  coefc[COEF_SIZE][NPOINT_APER/4];  /*  c  coef.  for  each  section  */ 
float  coefa[ COEF_SIZE]  [MAXDEG+ 1  ] ;  /*  a  coef.  for  each  c  */ 

long  deg[COEF_SIZE];  /*  degree  for  each  coefc  */ 

long  mat_c_stride,findex_c_stride,coef_c_stride; 

float  tstart[NPOINT_APER];  /*  time  from  trigger  to  start  data  acq.  at  */ 

/*  each  position  */ 

A.6.  Matrix  Pseudo  Inverse  Coefficients  (initmat.c) 

#include  <stdio.h> 

#include  <math.h> 

#include  "coefgen.const" 

#include  "coef gen. var” 

initmatO 

{ 

int  i,j; 
float  a,b; 

a=(  (float)  BOX_S  IZE_AZIMUTH-4.)/4. ; 
b=((float)BOX_SIZE__RADIAL-3.)/4.; 

mat[0]  [0]  =93 ./ 1 75 , ; 
mat[0]  [  1  ]=27 ./1 7 5 . ; 
mat[0][2]=(-9.)/175.; 
mat[0][3]=(-3.)/35.; 
mat[0][4]=9./175.; 
mat[0][5]=62./175.; 
mat[0]  [6] = 1 8./1 75.; 
mat[0][7]=(-6.)/175.; 
mat[0][81=(-2.)/35.; 
mat[0][9]=6./175.; 
mat[0]  [  1 0] =3 1 ./ 1 7  5 
mat[0][l  1]=9./175.; 
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mat[0][12]=(-3.)/175.; 

mat[0][13]=(-l.)/35.; 

mat[0][14]=3./175.; 

mat[0][15]=0.; 

mat[0][16]=0.; 

mat[0][17]=0.; 

mat[0][18]=0.; 

mat[0][19]=0.; 

mat[0][20]=(-31.)/175.; 

mat[0][21]=(-9.)/175.; 

mat[0]  [22]  -3./175.; 

mat[0][23]=l./35.; 

mat[0][24]-(-3.)/175.; 

tnat[l][0]==(-8L)/(175.*a); 

mat[l][l]=39./(350.*a); 

mat[  1  ]  [2] = 1 2  ./(3  5 .  *a) ; 

mat[l][3]=81./(350.*a); 

mat[l][4]=(-39.)/(175.*a); 

mat[l][5]=(-54.)/(175.*a); 

mat[l][6]=13./(175.*a); 

mat[l][7]=8./(35.*a); 

mat[l][8]=27./(175.*a); 

mat  [  1  ]  [9] =(-26. )/( 1 75 .  *a) ; 

mat[l][10]=(-27.)/(175.*a); 

mat[l][ll]=13./(350.*a); 

mat[l][12]=4./(35.*a); 

mat[l][13]=27./(350.*a); 

matr  1  ]  [  1 4]  — (~  1 3 .  )/(l  75.  *a); 

mat[l][15]=0.; 

mat(l][16]=0.; 

mat[l][17]=0.; 

mat[l][18]=0.; 

mat[l][19]=0.; 

mat[lj[20]=27./(175.*a); 

mat[l][21]=(-13.)/(350.*a); 

mat[l)[22]=(-4.)/(35*a); 

mat[l][23]=(-27.)/(350.*a); 

mat[l][24]=13./(175.*a); 

mat[2]  [0] =3  ./(3  5 .  *po  w  (a  ,2. ) ) ; 
mat[2][l]=(-3.)/(70*pow(a^.)); 

mat[2][2]=(-3.)/(35.*pow(a,2.)); 

mat[2][3]=(-3.)/(70.*pow(a,2.)); 

matr2][4]=3./(35.*pow(a,2.)); 

mat[2][5]=2./(35.*pow(a,2.)); 
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rnat[2][6J=(-l.)/(35.*pow(a,2.)); 
mat[2][7]=(-2.)/(35.*pow(a,2.)); 
mat[2][8]=(- 1  ,)/(35.  *pow(a,2.)); 
mat[2]  [9]  =27(35 .  *pow  (a,2 .)) ; 
cnatf2][10]=l./(35*pow(a,2.)); 
mat[2]  [  1 1  ]  =(- 1 .  )/(70.  *pow(a,2.)); 
mat[2][  1 2]=(- 1  .)/(35.  *pow(a,2.)); 
mat[2]  [  1 3]  =(- 1  .)/(70.  *pow(a,2.)) ; 
mat[2]  [  1 4]=1 7(35 .  *pow(a,2.)); 
mat[2][15]=0.; 
mat[2][16]=0.; 
mat[2][17]=0.; 
mat[2][18]=0.; 
mat[2][19]=0.; 

mat[2][20]=(-l.)/(35*pow(a,2.)); 
mat[2][21]=l./(70.*pow(a,2.)); 
mat[2][22]=l./(35.*pow(a,2.)); 
mat[2][23]=l  7(70.*pow(a,2.)); 
mat[2][24]=(-l.)/(35.*pow(a,2.)); 

mat[3][0]=(-31.)/(175.*b); 

mat[3][l]=(-9.)/(175.*b); 

mat[3][2]=37(175.*b); 

rnat[3][3]==l./(35.*b); 

mat[3][4]=(-3.)/(175.*b); 

mat[3][5]=(-31.)/(350.*b); 

mat[3][6]=(-9.)/(350*b); 

mat[3][7]=3./(350.*b); 

mat[3][8]=l./(70.*b); 

mat[3j[9]=(-3.)/(350.*b); 

mat[3][10]=0.; 

mat[3][ll]=0.; 

mat[3][12]=0.; 

mat[3][13]=0.; 

mat[3][14]=0.; 

mat[3][15]=317(350.*b); 

matf3][16]=9./(350.*b); 

mat[3][17]=(-3.)/(350.*b); 

mat[3]  [  1 8]  =(- 1 .  )/(70.  *b) ; 

mat[3][19]=3./(350.*b); 

mat[3][20]=31./(175.*b); 

mat[3][21]=9  ./(175. *b); 

mat[3][22]=(-3.)/(175*b); 

mat[3][23]=(-l.)/(35.*b); 

mat[3][24]=3./(175.*b); 


mat[4][0]=27./(175.*a*b); 

mat[4][l]=(-13.)/(350*a*b); 

mat[4][2]=(-4.)/(35.*a*b); 

mat[4][3]=(-27.)/(350.*a*b); 

matf-4}  [4]=1 3  V(  1 75.  *a*b); 

mar[4][5]=27./(350.*a*b); 

mat[4]  [6]=(- 1 3.)/(700  *a*b); 

mat[4][7M-2.)/(35.*a*b); 

mat[4][8]=(-27.)/(700.*a*b); 

inat[4][9]=13./(350.*a*b); 

mat[4]fl0]=0.; 

mat[4][ll]=0.; 

mat[4][12]=0.; 

mat[4][13]=0.; 

mat[4][14]=0.; 

mat[4]  [  1 51=(-27.)/(350.*a*b); 
mat[4]  [  1 6] = 1 3 ./(7 00.  *  a*b) ; 
mat[4][17]=2./(35.*a*b); 
mat[4][l  8]=27./(700.*a*b); 
mat[4][19]=(-13.)/(350.*a*b); 
mat[4][20]=(-27.)/(175.*a*b); 
mat[4]  [2 1  ]= 1 3 7(350.  *  a*b) ; 
mat[4][22]=4./(35.*a*b); 
mat[4][23]=27./(350.*a*b); 
mat[4]  [24]=(- 1 3.)/(175.*a*b); 

mat[5][0]=(-l.)/(35.*pow(a,2.)*b); 

mat[5][lM./(70.*pow(a,2.)*b); 

mat[5][2]=l./(35.*pow(a,2.)*b); 

mat[5  ]  [3] = 1  ./(70.  *po  w(a,2. )  *b); 

mat[5][4]=(~l.)/(35.*pow(a,2.)*b); 

mat[5]  [5]  =(- 1  .)/(70.*pow(a,2.)*b); 

mat[5][6]=l./(140.*pow(a,2.)*b); 

mat[5][7]=l./(70.*pow(a,2.)*b); 

mat[5][8]=l./(140.*pow(a,2.)*b); 

mat[5][9M-l.)/(70*pow(a,2.)*b); 

mat[5][10]=0.; 

xnat[5]fll]=0.; 

mat[5][12]=0.; 

mat[5][13}=0.; 

mat[5][14]=0.; 

mat[5]  [  1 5] =1  ./(70.  *pow(a,2.)*b); 
raat[5][  1 6]-(- l-)/(l  40.  *pow(a,2.)*b); 
mat[5][17]=(- 1  .)/(70.*pow(a,2.)*b); 
raat[5]  [  1 8] =(- 1 . )/( 1 40.  *pow(a,2. )  *  b) ; 
mat  [5]  [  1 9] = 1  ./(70.  *pow(a,2. )  *b) ; 


mat[5][20]=l./(35.*pow(a,2.)*b); 
mat[5][21  ]=(-!. )/(70.*pow(a,2.)*b); 
rnat[5][22]=(-l.)/(35.*pow(a,2.)*b); 
oiat[5][23]=(-l.)/(70.*pow(a,2.)*b); 
mat[5][24]=l./(35.*pow(a,2.)*b); 


) 


A.7.  Initialize  All  variables  (varinitc) 


J  ******  *  **  *  *  *  *  *  **  *  **  *  ***  *  ****  *  ******  ijp  *  *  *  *  **  *  ***  *  *  **  *  *  *  *  *  *  **  SH  sfc  *  *  *  *  *  sfc  *  *  *  *  *  :f:  *  j|c  * 


Subroutine:  coefgen_varinit 

Description  :  initialize  all  neccessary  variables 

sfc***************************************************************************  j 


#include  <stdio.h> 

#include  <math.h> 

#include  "coefgen. const” 

#include  "coefgen. var” 
coefgen_vaiinit() 

{ 

n_azimuth=NPOIbTT_AZIMUTH;  /*  #  of  points  in  azi.  direction  */ 

n_radial=NPOINT_RADIAL;  /*  #  of  points  in  rad.  direction  */ 

pix_size=n_azimuth*n _radial;  /*  Pixel  array  size  */ 

n_data=NPOINT_DATA;  /*  #  of  orig.  data  points  */ 

n_data_inter=NPOENTJDATA_INTER;  /*  #  of  interpolated,  data  points*/ 
n_data_inter_half=NPOINT_DATA_INTER/2;  /*  1/2  #  of  inter,  data  points  */ 
n  _aper=NPOINT_  APER;  /*  #  of  points  in  the  aperture  */ 

n_boxsample=NPOINT_BOXSAMPLE;  /*  #  of  samples  in  box  for  curve  fit  */ 
n_coef  =COEF_S  1ZE;  /*  #  of  coeffs  for  curve  fit  */ 

n  _section=SECTION . SIZE;  /*  #  of  sections  for  one  aperture*/ 


/*  Radar  Variables  */ 

fs=SAMPLING_RATE;  /*  Sampling  Frequency  */ 

ts=SAMPLING ^PERIOD;  /*  Sampling  period  */ 

c=3.e8;  /*  Speed  of  wave  */ 


/*  Geometry  variables  */ 

a_ofset=A_OFSET*PI/l  80.;  /*  Offset  angle  from  center  line  */ 

a_span=A_SPAN*PI/l80.;  /*  Spanning  angle  of  the  patch  */ 

d_jheta=a„span/((float)n_azimuth-  1 . );  /*  delta  theta  */ 

aper_length=(n_aper-l)*2.0*.0254;  /*  length  of  aperture  */ 

radar_height=RADAR_HEIGHT*  1 2.  *0.0254;  /*  height  of  radar  */ 

radar..  height2=radar_..height*radar  .height;  /*  square  the  height  */ 
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muncenter-RMINCENTER*12.*0.0254;  f*  min  range  from  center  position 

to  reference  point  */ 

d_rcenter=ts*c/(2.0*2.0);  /*  sampling  range  */ 

/*  project  back  to  2*n_data  samples  */ 
dr8=ts*c/(2. 0*8.0);  /*  interpolated  16K  buffer  */ 


/* 

theta„center=a_ofset-a_span/2.+((float)n_azimuth/2.-l.)*d_theta; 

c_theta_center=cos(theta_center); 

s_theta_center=sin(theta . center); 

reenter jrefcnmincenter+d_rcenter*((float)n_radial/2.-l.); 

*/ 

rcenter_ref=RCENTERJREF*  1 2. *.0254;  /*  range  of  ref.  point  from  center  pos.  */ 
theta_center_ref=THETA_CENTER„REF*PI/180.;  /*  angle  of  ref.  point  from  center 

position  */ 

c_theta„center_ref=cos(theta_center_ref); 

s_theta_center  jef=sin(theta_center_ref); 


/*  SSL  Variables  */ 

f_zero=0.;  /*  floating  point  zero  */ 

stride=l;  /*  stride  for  ssl  */ 

f_one=1.0;  /*  floating  point  1  */ 

i_one=l;  /*  integer  1  */ 

mat_c_stride=n„boxs ample;  /*  column  stride  for  mat 

findex  c . stride— 1 ;  /*  column  stride  for  findex 

coef_c_stride=l;  /*  column  stride  for  coeff 

deg[0]=3;  /*  degree  for  first  coefc  */ 

deg[l]=3; 

deg[2]=3; 

deg[3]=3; 

deg[4]=2; 

deg[5)=2; 


*/ 


*/ 

*/ 


/*  xgints„(pixel,&f_zero,&pix_size,&stride);  *//*  clear  pixel[][]  */ 

} 

A.8.  Find  Least  Squares  Polynomial  Fit  (poly  fitc) 

^***** ********************************************************** ********** 
Subroutine:  poly_fit.c 
Description: 

Perform  data  fit  to  polynomial 
Input:  x[n]  input  array 
y[n]  input  array 

n  number  of  elements 

deg  degree  of  poly. 
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output  coeff[deg+l]  coefficients  of  poly 

*  **************  *********  ifc************  Jlc****************  Jit****:):******  *:}:*  ***•!(:  :jy 

#include  <stdio.h> 

#include  <malloc.h> 

#include  <math.h> 


void  polyjfit(x,y,n,deg,coeff) 

float  x[],y[],coeff[]; 
long  n,deg; 

{ 

double  *dx,*dy,*dcoeff; 

double  *dX,*dXtrans  *dA,*dB  *dC; 

long  m,ij,one; 

one=l; 


m=deg+l; 

dx=malloc(n*sizeof(double)); 

if(dx==NULL) 

{ 

printf("Memorv  allocation  Error  !!!\n"); 
exit(l); 

} 

dy=malloc(n*sizeof(double)); 
if(dy— NULL) 

{ 

printf(”Memory  allocation  Error  !!!\n’'); 
exit(l); 

} 

dcoeff=malloc(m*sizeof(doubIe)); 

if(dcoeff==NULL) 

{ 

printff'Mernory  allocation  Euor  !!!\n"); 
exit(l); 

} 

dX=malloc(n*m*sizeof(doubIe)); 
if (dX==NULL) 

{ 

printf("Memory  allocation  Error  !!!\n"); 
exit(l); 

} 

dXtrans=malloc(n*m*sizeof(double)); 

if(dX==NULL) 

{ 

printf("Memory  allocation  Error  !!!\n"); 
exit(l); 

} 

dA=maUoc(m*m*sizeof(double)); 

if(dA==NULL) 
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{ 

printf("Memory  allocation  Error  !!!\n"); 
exit(l); 

) 

dB=malloc(m*m*sizeof(double)); 

if(dB==NULL) 

{ 

printfC'Memory  allocation  Error  ! !  !\n 
exit(l); 

) 

dC=malloc(m*sizeof(double)); 

if(dx==NULL) 

{ 

printfC'Memory  allocation  Error  H!\n"); 
exit(l); 

} 

for(i=0;i<n;i++) 

{ 

dx[i]=(double)x[i]; 

dy[i]=(double)y[i]; 

} 

for(i=0;i<n;i++) 

{ 

dX[i*m]=1.0; 

) 

for(j=l;j<m;j++) 

{ 

for(i=0;i<n;i++) 

{ 

dX  [j  +i  *m]  =pow  (d  x  [i] ,  (dou  ble)j) ; 

) 

} 

mxtrans(dX,n,m,dXtrans); 

mxmuld(dXtrans,&n,&one,dX,&m,&one,dA,&m,&one,&m,(&m,&n); 

mat_inverse(dA,dB,m); 

mxmuld(dXtrans,&n,&one,dy,&one,&one,dC,&one,&one,&ni,&one,&n); 

mxmuld(dB,&m,&one,dC(&one,&one,dcoeff,&one,&one,&m,&one,&m); 

for(i=0;i<m;i++) 

{ 

coefffi] =(float)dcoeff[i] ; 

} 

free(dx); 

free(dy); 

free(dA); 

free(dB); 
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free(dC); 

free(dX); 

free(dXtrans); 

free(dcoeff); 

} 


A.9.  Matrix  Inverter  (inverse.c) 


#include  <stdio.h> 

mat_inverse(source_mat,dest_mat,size_mat) 
double  source_mat[]; 
double  dest_mat[]; 
long  size_mat; 

{ 

long  n,i,j,k,ki,count; 
double  b,bl; 
double  err=0.0001; 
double  a[100][100]; 
n=size_mat; 
count=0; 

for(i=0;  i<size_mat;  i++) 

{ 

for(j=0;j<size_mat;j++) 

{ 

a[i+ 1  ]  [j+ 1  ]=source_mat[count++] ; 

} 

} 

for(i=l;i<=n;i++) 

{ 

for(j=n+ 1  ;j<=2*n;j++) 

{ 

if((j*n-i)==0) 

{  a[i][j]=1.0;  } 
else 

{  a[i]|j]=0.;  } 

} 

} 

for(k=  1  ;k<=n- 1  ;k++) 

{ 

b=a[k][k]; 
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ki=k; 

for(i=k+l;i<=n;i++) 

{ 

if((abs(b)-abs(a[i]  [k]))<0) 

{ 

b=a[i][k]; 

ki=i; 

} 

} 

if((abs(b)-err)<0) 

{ 

printf("Error  matrix  inverse  ,  matrix  is  singular\n"); 
exit(l); 

} 

if((ki-k)!=0) 

{ 

for(j=k;j  <=2*n;j++) 

{ 

bl=a[k][j]; 

a[k][j]=a[ki][j]; 

a[ki][j]=bl; 

} 

} 


for(j=k+l;j<=2*n;j++) 

{ 

a[k][j]=a[k][j]/b; 

} 

for(i=k+l;i<=n;i++) 

{ 

for(j=k+ 1  ;j<=2*n;j++) 

{ 

a[i]Q]=a[i][j]-a[i][k]*a[k][j]; 

} 

} 

} 

for(j=n+l;j<=2*n;j++) 

{ 

a[n]  [j]=af  n](j]/a[n]  [nj; 


for(k~n- 1  ;k>=  1  ;k=k- 1 ) 

{ 
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for(j=n+l;j<=2*n;j++) 

{ 

for(i=k+ 1  ;i<=n;i ++) 

{ 

a[k]  [j  a[k](j]-a[k]  [i]  *a[i]  [j]; 

} 

} 

} 

count=0; 

for(i=0;i<size_mat;i++) 

{ 

for(j=0;j<size_mat;j ++) 

{ 

dest_mat[count++]=a[i+ 1  ]  [size_mat+j+l]; 

} 

} 
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APPENDIX  B 


B.l.  Main  Program  For  SUN  Computer(focus.c) 


program:  focus.c 
description: 

this  programs  reads  radar  data  from  file,  interpolates 
data  using  convolution  in  freq.  domain,  and  projects 
data  back  to  the  polar  grid  on  the  ground.  The  indeces 
used  in  back  projection  are  computed  from  the  file 
which  contains  coefficients  for  a  particular  geometry 


#include  <stdio.h> 

#include  "focus.h” 

#include  "focus. var" 

float  da  ta[NPOINT_D  AT  A] ;  /*  data  buffer  for  host  */ 
struct  cstype  *csfNO_SUPERCARD]; 

main  () 

{ 

float  tstart_err[NPOINT_APERj;  /*  tstart-tstart_trunc  for  each  position  */ 
FILE  *inputfile;  /*  input  file  contains  radar  data  */ 

FILE  *outputfile;  /*  output  file  contains  focused  image  */ 

FILE  *filterfile;  /*  input  file  contains  filter  coefficients  */ 

FILE  *coeffile;  /*  input  file  contains  back  projection  */ 

/*  coefficients  */ 

FILE  *tstart_err_file;  /*  input  file  contains  tstart-tstart  trunc  */ 

FILE  *junkfile; 
char  inputname[80]; 
char  outputname[80j; 
char  filtemame[80]; 
char  coefname[80j; 
char  tstart_err_name[80]; 
char  temp_string[20]; 

long  sec tion ,pos ition, box,c oefc_order,coefa  order; 
longij;  /*  just  working  variable  */ 
float  ^position; 

long  data_ready=l;  /*  maibox  to  indicate  data  ready  host  —  >  SC  */ 
long  data_consumed=2;  /*  mailbox  to  indicate  data  consumed  host  <—  SC  */ 
long  msg;  /*  message  read  from  mailbox  */ 
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long  one=l;  /*  nonzero  message  to  put  in  mailbox  */ 
long  cacheset=l;  /*  cacheset=0  turn  off  cache  */ 
long  numpar=l; 

float  pixel_zero[BOX_SIZE_AZIMUTH|  [NPOINTRADIAL] ; 


!*#*  open  supercards  *********************************************/ 
for  (i=0;i<NO„SUPERCARD;i++) 

{ 

cs[i]=(struct  cstype  *)xlubgn_(0,&cacheset,"sc.lo"); 

} 


yfjjfc c jfe  SUpCTCHXd  V331iiblcs 

focus_varinit(); 


/***  initialize  supercards’  id  ***********************************/ 
for(i=0;i<NO_SUPERCARD;i++) 

{ 

cs  [i]->supercard_id=i; 

} 


Enter  Input  ’t'************************************************/ 
printf("Enter  input  file  name  ( .raw  )==>  "); 
scanf("%s’\inputname); 
printf("Enter  output  file  name  ( .focus)===> 
scanf(”%s",outputname); 

/* 


printf("Enter  back  projection  coefficient  file  name  (coefa.dat)===>  "); 
scanf("%s",coefname); 

printf("Enter  start  time  difference  file  name  (delaytime_diff.dat)~=>  "); 
scanf("%s",tstart_err„name); 

*j 


strcpy(coefname,''coefa.dat"); 
strcpy(tstart_err_name, "delaytime_diff.dat"); 


printf("Do  you  want  hamming  weight  accross  the  aperture  (y  or  n)  ===>  "); 
scanf("%s",temp_string); 


for(j=0;j<NO__SUPERCARD;j++) 

{ 

csfj]  -  >ham_flag=0; 

if(temp . string[0]  ==’y’ ) 

{ 

cs[j]->ham_flag=l; 

} 
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strcpy(filtemame,  "filter.dat"); 
if((inputfile=fopen(inputname,"rb”))-=NULL) 

{ 

printf(" Error  open  input  file\n"); 
exit(l); 

} 

if((outputfile=fopen(outputname,,'wb"))==NULL) 

f 

printf("Error  open  output  file\n"); 
exit(l); 

} 

if((filterfile=fopen(filternarne,"r"))==NULL) 

{ 

printf("Error  open  filter  file\n"); 
exit(l); 

} 

if((coeffile=fopen(coefname,"r"))==NULL) 

{ 

printf("Error  open  back  projection  coefficient  fileXn"); 
exit(l); 

} 

if((tstart_err_file=fopen(tstart_err_name,"r"))==NULL) 

{ 

printf("Error  open  tstart  difference  fileNn"); 
exit(l); 

} 

/***  filter  coefficients  **********************************/ 
/***  into  1st  supercard  and  copy  to  others  supercards  **********/ 
printf("»>  Reading  filter  coefficientsW); 
for(i=0;i<FILTER_LENGTH;i++) 

{ 

fscanf(filterfile,"%f',&cs[OJ->filter_coe^iJ); 

) 

for(j= 1  ;j<NO_SUPERCARD;j++) 

{ 

for(i=0;i<FILTER_LENGTH;i++) 

{ 

cs[j]->filter_coef[i]=cs[0]->filter_coef[i]; 

} 

} 
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/***  Read  in  tstart_enr  for  all  positions  to  first  supercard  *****/ 

/***  and  then  copy  to  other  supercards  ***************************/ 
printf("»>  Reading  tstart__err  data  \n"); 
for(i=0;i<NPOINT_APER;i++) 

{ 

fscanf(tstart . err_file,"%d  %f',&j,&(cs[0]->tstart_err[i])); 

} 

for(j= 1  ;j<NO  SUPERCARD;j++) 

{ 

for(i=0;i<NPOINT_APER;i++) 

{ 

cs  [j]->tstart„err[i]  =cs[0]->tstart_err[i] ; 

} 

} 


/***  Read  in  coefficients  a  for  back  projection  ******************/ 

/***  to  the  first  supercard  and  then  copy  to  other  supercards  ****/ 
printf("»>  Reading  backprojection  coefficients  \n”); 
for(section=0;section<SECnON_SIZE;section++) 

{ 

for(box=0;box<NBOX_RAD3AL*NBOX„AZIMUTH:box++) 

( 

for(coefc„order=0;coefc_order<COEF_SIZE;coefc__order++) 

{ 

for(coefa__order=0;coefa_order<(cs[03->deg[coefc_order]+l);coefa_order++) 

{ 

fscanf(coeffile,"%f',&(cs[0]->coefa[section][box][coefc_order][coefa__order])); 
for(j=l  ;j<NO_SUPERCARD;j++) 

{ 

cs  [j]  ->coefa[section]  [box]  [coefc_order]  [coefa_order] = 
cs[0]  ->coefa[section][box]  [coefc_order]  [coefa_order] ; 

} 

}  /*  coefa_order  */ 

}  /*  coefc__order  */ 

}  /*  box  */ 

}  /*  section  */ 


/***  supercards  ********************************************/ 
printf("»>  Starting  supercard  programs  \n"); 
for(j =0;j  <N  O _S  UPERC ARD ; j ++) 

{ 

xrcaU_("sc"J&numpar,&csp]->dummy); 

} 


/***  For  every  position,  the  host  computer  reads  data  for  that  ***/ 
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/***  position  and  copy  data  into  each  supercard  memory.  **********/ 

/***  There  are  two  mail  boxes  used  betwwen  host  and  each  *********/ 

/***  supercard  to  provide  synchronization  ************************/ 

for(section=0;section<SECTION„SIZE;section++) 

{ 

for(position=0;position<NPOINT_APER/SECTION_SIZE;position++) 

{ 

printf("»»>processing  section  %d  pos  %d  «««<\n",section,position); 

fread(data,sizeof(data),  1  ,inputfile); 

for(j=0;j<NO . SUPERCARD;j++) 

{ 

/***  Copy  data  in  each  supercard  memory  **************************/ 
for(i=0;i<NPOINT_DATA;i++) 

{ 

cs  fj  ]  -  >data. ..temp  [i]  =data[i] ; 

} 

/***  Send  mail  to  supercard  to  inform  data  is  ready  **************/ 
xlnwxmt_(cs  (Jj  ,&data_ready  ,&one); 

} 

/***  Wait  until  data  is  consumed  by  supercards  *******************/ 
for(j =0;j  <NO„S  UPERC  ARD ; j ++) 

{ 

xlwtrec_(csfj],&data... consumed, &msg); 

} 

}  /*  position  */ 

}  /*  section  */ 


/***  Check  for  supercards  to  finish  all  processing  ******************/ 

for(j=0;j<NO . SUPERCARD;j++) 

{ 

xldone„(cs[j]); 

} 

!***  Save  resulting  radar  image  *************************************/ 
for(j=0;  j  <BOX_B  ASE_START ;  j  ++) 

{ 

printf("»>Prepad  zero  to  output  file\n"); 
fwrite(pixel_zero,sizeof(pixel_zero),l,outputfile); 

} 

printf(">»Save  image  to  output  file\n"); 
for(j=0;j<NO_SUPERCARD;j++) 
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{ 

fwrite(cs|j']->pixel, 

sizeof(float)*BOX_SIZE_AZIMUTH*NPOINT_RADIAL*NO_KBOX_PROCESS/NO_SUPE 

RCARD, 

l,outputfile); 

} 

for  (j  =B  OX_B A SE_ST ART+NO_KBOX  PROCESS;j<NBOX_AZIMUTH;j++) 

{ 

printf("»>Postpad  zero  to  output  file  \n"); 
fwrite(pixel_zero,sizeof(pixel  zero),  1  ,outputfile); 

} 


fclose(inputfile); 

fclose(outputfile); 


/***  Closing  all  supercards  ****************************************/ 
for(j=NO_SUPERCARD- 1  ;j>=0;j~) 

{ 

xlclos_(cs(j]); 

} 


}  /*  end  main  */ 


B.2.  Main  Program  For  Multiple  Array  Processors  (sc.c) 


y* ****************************************************** ************ 


Program:  sc.c 


Description:  This  is  the  main  program  to  be  executed  by  each  of 
the  CSPI  i860  array  processor. 


*****************************************************^ 


#include  "focus.h'1 
#include  "focus.var" 


struct  cstype  *cs; 


void  sc(dummy) 
long  *dummy; 

{  ' 

void  focus_varinit(); 
void  filter_init(); 
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void  haminit(); 

void  hamwt(); 

void  eraseQ; 

void  inter(); 

void  fix„data_pointer(); 

void  poly(); 

void  bp(); 

void  xwtrec_(); 

void  xnwxmt_(); 

void  xvmov_(); 

void  xvclr_(); 

float  *data_pointer;  /*  pointer  to  actual  data  */ 

long  section,position,box,coefc_order,kbox,ibox; 

long  data . ready=l;  /*  maibox  to  indicate  data  ready  host  — >  SC  */ 

long  data_consumed=2;  /*  mailbox  to  indicate  data  consumed  host  <—  SC  */ 

long  msg;  /*  message  read  from  mailbox  */ 

long  one=l;  /*  nonzero  message  to  put  in  mailbox  */ 

float  f_position; 

long  boxbase; 

initialize  supercard  pointer  ***********************************/ 
cs=0; 


initialize  image  with  zeros  ************************************/ 
xvclr_(cs->pixel,&cs->pix_size,&cs->i„one); 

-*fc  3^C  2^c  5^C  }Jc  jjc  j|c ^  ijc ijc  jjc  jfv  3^C 5^ 

fllter_init(); 

/***  initialize  hamming  weight  coefficients  *************************/ 
haminit(cs->ham_coef,NPOINT_APER); 

gtart  f]  ^  lniH^C  5k  5jc  ^  iff  5i»  jf*  ?k  »f»  jft  ?f<  4s ‘f1*  %  ^  ^  ^  'f' ^  'i' 'J*  ^ 

boxbase=BOX_BASE„START+cs->supercard_ids!tNO_KBOX_PROCESS/NO_SUPERCARD; 

for(section=0;section<SECTION„SIZE;section++) 

{ 

for(position=0;position<NPOINT_APER/SECTION„SIZE;position++) 

{ 

fox*  new  data,  fxom  host  computer 
x  wtrec_(&data_ready,&msg) ; 

/***  Copy  data  to  working  buffer  ***********************************/ 
xvmov_(cs->data,cs->data_temp,&cs->n_data,&cs->Lone,&cs->i_one); 
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/***  Signal  host  computer  that  data  are  consumed  *******************/ 
xnwxmt_(&data_consumed,&one); 

f_position=(float)position; 

/***  Hamming  weight  data  if  needed  *********************************/ 
if(cs->ham_flag==l) 

{ 

hamwt(cs->data,cs->n_data,section*NPOINT_APER/SECTTON__SIZE+position, 

cs->ham_coef); 

} 

/***  zeros  last  portion  of  data  for  circular  convolution  ***********/ 
erase  (c  s-  >data,cs->n_data,  32) ; 

Interpolate  St**^!^******************!^*****)!:!):*^!******^;**^!*^:*^ 

inter(cs->data,cs->n_data,cs->data . inter, 

cs->n_data . inter,cs->filter_fft); 

fix_data_pointer{&data_pointer,&cs->data_inter[0], 

cs->tstart_err[NPOINT_APER/SECTION . S  IZE*section+position] , 

cs->ts,8); 

for(kbox=boxbase;kbox<(boxbase+NO_KBOX_PROCESS/NO_SUPERCARD);kbox++) 

{ 

for(ibox=IBOX_START;ibox<(IBOX_START+NO_IBOX_PROCESS);ibox++) 

{ 

box=kbox  *NB  OX_RAD  IAL+ibox; 

/***  Generate  coefficients  for  back  projection  ********************/ 
for(coefc_order=0;coefc_order<COEF_SIZE;coefc_order++) 

{ 

poly(&f_position,&cs->coefc[coefc_order],  1 , 

&cs->coefa[section]  [box]  [coefc_order]  [0] , 
cs->deg[coefc_order]); 

}  /*  coefc_order  */ 

I***  Perform  {jack  projection  fc*************************************/ 
bp(cs->pixel,data_pointer,cs->n_data_Jnter,cs->coefc, 
boxbase  ,kbox  ,i  box, 

cs->k_pix_size,cs->i_pix„size,cs~>i_box_size); 

}  /*  i_box  */ 

}  /*  k„box  */ 

}  /*  position  */ 
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}  /*  section  */ 

}  /*  end  supercard  program  */ 

B.3.  Routines  Used  to  Interpolate  (inter.c) 


#include  <math.h> 
#include  "focus.h" 
#include  "focus.var" 


y****  ******  **4:  **************  *********************  *******  *********  ******  ****** 


Subroutine:  haminitx 


Description: 

this  subroutine  initialize  array  of  hamming  weighting  coeffs 
accross  the  aperture. 

Input:  float  ham_coef[NPOINT_APER] 

long  npoint  Number  of  points  for  hamming  coeffs 
Output:  float  ham_coef[NPOINT_APER]  is  initialized 


haminit(ham_coef,npoint) 
float  *ham_coef; 
long  npoint; 

{ 


long  i; 

for(i=0;  icnpoint;  i++) 

{ 

ham_coef[i]=.54-.46*cos(2.*3.1415927*(float)i/(float)npoint); 

} 


} 


y**************** ************************************************************ 


Subroutine:  hamwt.c 


Description: 

this  subroutine  takes  the  radar  signal  at  a  position  and 
multiplies  the  entire  signal  array  with  the  hamming  coef  at 
that  position. 

Input: 

float  data[NPOINT_DATA] 
long  npoint_data 
long  position 
float  ham_coef[]; 


Output: 

float  data[NPOINTJDATA] 


72 


£;£*****  ***  ********  ***********  ****  **  **********  *******************************/ 


hamwt(data,npoint„data,position,ham_coef) 
float  *data,*ham_coef; 
long  npoint_data, position; 

{ 

long  i; 

f or(i=0;  i<npoint_data;i++) 

{ 

data[i]  =data[i]  *h  am__coef[position] ; 

} 

} 


j  ***  *  *****  *  *******  *  *  *  **************  *  ************  *  ***  *  *  **  *  *  ********  *  ***  *  *****  * 


Subroutine:  inter. c 
Description: 

This  routine  performs  FIR  interpolation  by  using  convolution 
in  the  frequency  domain.  The  input  buffer  contains  2K  of  data 
and  is  interpolated  into  16K  of  data.  The  FIR  filter  has 
breakpoints  at  .95  Ghz  and  1.05  Ghz 
#  of  taps  =  255 

name  of  filter  coeff.  file  :  filter.dat 


Input:  data[  NPOINT_D  AT  A] 

complex  filter_fft[NPOINT_DATA_INTER/2J 
float  nyquist_filter 

Output:  data_inter[NPOINT_DATA,INTER] 


inter() 

{ 

void  xvclr_(); 
void  xfrf_(); 
void  xcvmls_(); 
void  xfri_(); 

extern  struct  cstype  *cs; 
int  i; 

/*  clear  data_inter  buffer  */ 

xvclr_(cs->data_inter,&cs->n_data_inter,&cs->i_one); 


/*  interleave  data  into  data-inter  buffer  */ 
for(i=0;i<NPOINT_DATA;i++) 

{ 
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cs->data_inter[i*83=cs->data[i]; 

} 

/*  FFT  of  interleaved  data  */ 

xfrf_(cs->data_inter_fft,&cs->nyquist_data,cs->data_inter, 

&cs->n_data_inter_half); 

/*  FFT  of  interpolated  data  (multiply  with  filter  in  freq.  domain)*/ 
xcvmls_(cs->data_inter_fft,&cs->f_one,cs->dataJnter_fft, 
cs->filter_fft,&cs->n_data_inter_half); 
cs->nyquist_data=cs->nyquist_data*cs->nyquist_fllter; 

/*  inverse  FFT  to  get  interpolated  data  */ 
xMJcs->dataJmer,&cs->nyquist_data,cs->data_inter_fft, 
cs->data_inter_fft,&cs->n_data_inter_half); 


} 


^  ^  ^  -f-  ^  ^  ^1*  ^  ^  ^  ^  ^  *^*  ij^  ^  ^  *|r  ^  i|^  <j|*  ^  4^ 


Subroutine:  fix_data_pointer 

Description:  This  subroutine  fixes  the  pointer  to  the  start  of 
data_inter  buffer.  The  pointer  needed  to  be  adjusted 
due  to  the  following  reasons: 

(a)  The  linear  FIR  interpolation  filter  has  phase 
shift  and  the  actual  data  point  starts  at  bin  127 
of  the  data_inter  buffer  ( FIR  has  255  coeffs  ) 

(b)  The  actual  time  to  start  data  acquisition  has  to 
be  rounded  of  to  1  nsec  resolution  for  the  scope 
DSA602.  However  the  backprojection  algorithm  uses 
the  exact  time  since  the  indeces  have  to  be  smooth 
for  curve  fit. 


Input: 

float  *data_inter  (  address  of  data_inter  buffer ) 
float  tstart_err  ( tstart-tstart_trunc  ) 
float  ts  ( radar  sampling  period  ) 
long  inter_factor  (interpolation  factor  ) 
Output: 

float  *data_pointer  ( actual  start  pointer  for  data  ) 


********************************************************:#:  *****************  j 
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flx_data_pointer(pointer,dataJnter,tstart_err,ts,inter_factor) 

float  **pointer; 

float  *data_jnter; 

float  tstart_err; 

float  ts; 

long  inter_factor; 

{ 

long  i; 

double  float l,float2; 

unsigned  long  index„ofset; 

float  l=modf(tstart_err*inter_factor/ts,&float2); 

if  (float  1  >==.  5 )  float2=float2+l.; 

index_ofset=(unsigned  long)127+(unsigned  long)float2; 

*pointer=data_Jnter+index_ofset; 

/*  zeros  fill  first  128  points  of  data_inter  since  phase  shift  from  linear  filter  */ 
for(i=0;i<  1 28;i++) 

{ 

data_inter[i]=0.; 

) 

/*  zeros  fill  the  last  128  points  of  data_inter  plus  8  points  for  tstart_err  */ 
for(i=0;i<136;i++) 

{ 

data_inter[NPOINT_DATA_INTER+i]=0.; 

} 

/*  first  point  and  last  point  of  actual  data  array  are  zeros  for  backprojection  */ 
**pointer=0.; 

*(*pointer+NPOINT_DATA_INTER- 1  )=0.; 

} 


^* ************** ********************************************************** 


Subroutine:  erase 

Description:  Zero  out  the  last  nzero  points  of  input  buffer 
************* **************************************************************/ 


erase  (b  uffer ,  ndata,  nzero) 
float  ^buffer; 
long  ndata, nzero; 

{ 
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long  i; 

for(i=ndata-nzero;i<ndata;i++) 

{ 

buffer[i]=0.; 

} 


} 


B.4.  Initialize  Interpolation  Filter  (filtinitc) 


y**********  *  ******  **  **  **********  sfc  =)=  *  ***************************  ****  ********* 


This  subroutine  performs  the  following  : 

1-  Read  filter  coefficients  from  file  to  filter  buffer 

2-  Zerro  pad  the  filter  buffer  to  the  size  of  interpolated  data 
buffer  size 

3-  Perform  Real  Forward  FFT  of  filter  buffer  to  prepare  for 
interpolation 


^include  <stdio.h> 
#include  "focus.h" 
#include  "focus.var" 


extern  struct  cstype  *cs; 

void  filter_init() 

{ 

void  xfrf_(); 
long  i; 

/*  zero  pad  filter  coeff  */ 

for(i=FILTER_LENGTH;i<NPOINT_DATA_INTER;i++) 

{ 

cs->filter_coef[i]=0.; 

} 

/*  take  fft  of  filter  coef  for  each  supercard  */ 
for(i=0;i<NO_SUPERCARD;i++) 

{ 

xfrf_(cs->filter_fft,&cs->nyquist_filter, 

cs->filter_coef,&cs->n_data__inter_half); 

} 


} 

B.5.  Main  Back  Projection  Focusing  Routing(bp.c) 

#include  <math.h> 
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y*:*  ******  ******  ********************************  ***  ******************  ********* 

Subroutine:  poly.c 
DEscription: 

Calculates  values  of  polynomials 


Input:  x[] 

n 

coeffl'j 

deg 

Output:  y[] 


input  array 
number  of  elements 
coefficients  of  poly 
degree  of  poly 
output  array 


t ****************************************************************** *********^r 


poly(x,y,n,coeff,deg) 
float  *x,*y,*coeff; 
long  n,deg; 

{ 

long  i  j; 

for(i=0;i<n;i++) 

{ 

y[i]=coeff[0]; 

for(j=l;j<(deg+l);j++) 

{ 

y[i]=y[i]+coeffO]*pow(x[i]((float)j); 

} 

} 


y* ******************************************************************* ******* 


Subroutine:  poly2.c 

Description:  calculate  y=c0+cl*x+c2*xA2 
where  x=[0,l,2,...,n-l] 


Input: 

float  c[3]  coefficients 
long  n  (max  :  1000)  number  of  points 
Output: 

float  y[]  output  vector 


#define  MAX.ELEMENT  1000 


void  poly2(y,n,c) 
float  y[],c0; 
long  n; 

{ 

void  xvrmp_(); 
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void  xdintg_(); 
float  temp; 

float  y  1  [MAX„ELEMENT] ; 
float  yl_init,yl_inc; 
yl_init=c[l]-c[2]; 
yl_inc=2.*c[2]; 

x  vrmp_(y  1  ,&temp  ,&y  1  _init,  &y  l„inc ,  &n) ; 
yl[0]=0.; 

xdintg_(y,&c[0],yl,&temp,&n); 

} 


j  3ft  jjc  3ft  ^fC  3ff  ifc  ijc  3^  3jt  3^  3|t  3|f  ifc  ?ft  %  3ft  #j(  jjf  3|C  3jC  3|t  5ft  ^  3fc  3fv  #|C  3ft  3ft  3ft  Sft  Sft  3ft  3jC  3ft  3ft  3ft  3|C  3f*  3ft  3ft  3$t  3ft  jjt  *ft  3ft  3ft  3ft  3ft  3ft  3ft  3ft  jjt  3ft  3ft  3jt  3f»  3ft  3$C  3ft  tjt  3ft  3ft  3ft  3ft  3ft  3ft  3ft 


Two-dimensional  back  projection  subroutine 

using  fast  algorithm  by  Nuttal  to  calculates  the  index  to  data 

array  for  each  pixel  on  the  ground  in  i_box,k_box  positioned  by  i,k 

k:  azimuth  (  2nd  order) 

i:  range  ( ist  order) 

The  equation  for  index  calculation  is 

index=  cO+cl*k+c2*kA2  +  (c3+c4*k+c5*kA2)*i 
=  dO  +  dl*i 

where  c(i)  (i=0,5)  is  a  set  of  coefficients  for  each  box  and 
each  position  in  the  aperture 
Input- 

float  *pixel  pointer  to  first  point  of  the  image 
2-dimensional  area 
float  *data  pointer  to  data  array 
Important:  data[0]=0.0 

datafdat_size- 1]=0.0 
long  dat_size  size  of  data  buffer 
float  *c  pointer  to  six  coefficients  for 
index  computation 

long  k_pix_size  size  of  the  patch  (pixels)  in  k  axis 
long  i_pix_size  size  of  the  patch  (pixels)  in  i  axis 
long  k_.  box  patch  number  in  k  axis 
long  i_box  patch  number  in  i  axis 
long  i_box_size  number  of  patches  in  i  axis 
Output: 

index  to  data  array  is  calculated  and  pixel  is  updated 


3fc  jfe  jft  3ft  3fc  3ft  jfc  3fC  3fC  jfc  3fc  3f33fc  3ffi  3fC  3ft  3fc  3fc  3fC  «fC  3fc  Jft  3ft  3fc  3ft  3fc  rft  3ft  3ft  3fc  3ft  tft  3f»  3fC  3fC  3fC  3fC  3fC  3fC  3ft  3fC  3fC  3ft  3fe  3ft  jft  jft  jft  3fc  jft  3fC  jfc  3fC  3ft  3fc  3ft  3ft  jft  jft  jft  3fc  3ft  3ft  3ft  jft  jft  3ft  jft  3ft  3ft  3fc  3ft  3ft  jft j 


#define  MAX_K_PIX_SIZE  1000 
#define  MAX_I_PIX_SIZE  1000 
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void  bp(pixel,data,dat_size,c,k_box_base,k_box,i_box,k_pix_size, 

Lpix_size,i_box_size) 

float  *pixel;  /*  array  of  pixels  of  the  whole  2-dimentional  area  */ 

float  *data;  /*  data  array  */ 

long  dat_size;  /*  size  of  data  array  */ 

float  *c;  /*  pointer  to  six  coefficients  for  index  computation  */ 

long  k_box_base;  /*  0  for  1st  SC,  i*NBOX_AZIMUTH/NO_SUPERCARD  for  ith  SC  */ 

long  k_box;  /*  patch  number  in  k  axis  */ 

long  i_box;  /*  patch  number  in  i  axis  */ 

long  k_pix_size;  /*  size  of  the  patch  (pixels)  in  k  axis  */ 

long  Lpix„size;  /*  size  of  the  patch  (pixels)  in  i  axis  */ 

long  i_box_size;  /*  number  of  patches  in  i  axis  */ 

{ 

void  xvclip_(); 
void  vclip_(); 
void  xvfx4_(); 
void  fix4_(); 
void  xvrmp_(); 
void  vindex_(); 
void  vramp_(); 
void  vgathr_(); 
void  vadd_(); 
void  vtabi„(); 

float  findexfM AX_I_PIX_S IZE];  /*  data  index  value  */ 

long  lindex[MAX_I_PIX_SIZE];  /*  data  index  value  (round  to  integer)  */ 

float  pix_temp  [MAX_I_PIX_S1ZE] ;  /*  temp  buffer  for  back  projection  */ 

float  dO[M AX_K_PIX_SIZE] ;  /*  d0=  cO+cl  *k+c2*kA2  */ 

float  dl  [M AX_K_PEX  J5IZE] ;  /*  dl  =  c3+c4*k+c5*kA2  */ 

float  *pix_index;  /*  pointer  to  current  pixel  */ 

long  pix_index_inc;  /*  each  time  k  is  incremented  ,  pixel  pointer  is  jumped  */ 

long  k;  /*  for  k  and  i  loops  control  */ 

float  maxindex;  /*  maximum  value  for  data  index  */ 

float  zero=0.0; 

float  one=  1.0; 

long  i_one=l; 

float  temp; 

maxindex=(float)(dat_size~  1 ) ; 

pix_index=pixel+(k„box-k_box_base)*k_pix_size*i_box_size* 

i_pix„size+i_box*i_pix_size; 

/*  index  of  first  point  of  the  patch  */ 
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/*  with  respect  to  pixel  */ 

pix„index_inc=i_pix_size*i_box_size;  /*  index  increment  along  k  axis  */ 

poly2(d0,k_pix_size,c);  /*  generate  array  of  dO  */ 

poly2(dl,k„pix_size,c+3);  /*  generate  array  of  dl  */ 

for(k=0;k<k_pix_size;k++) 

{ 

/*  generate  floating  point  index  vector  */ 
vramp_(&dO[k]  ,&d  1  [k]  ,fmdex,&i_one,&i_pix_size); 

/*  table  look  up  and  interpolate  */ 

/*  this  function  replaces  xvclip(),xvfx4(),vgathr()  */ 

/* 

vtabi_(findex,&i_one,&one,&zero,data,pix_temp,&i_one,&dat_size,&i_pix_size); 

*/ 

/*  clip  index  */ 

vclip_(findex,&i_one,&zero,&maxindex,findex,&i . one,&i . pixsize); 

/*  convert  to  integer  index  */ 
fix4_(findex,&i_one,lindex,&i_one,&Lpix_size); 

/*  gather  data  into  temporary  buffer  */ 
vgathr_(data,lindex,  &i_one,pix_temp,&Lone,&i_pix_s  ize) ; 

/*  update  pixel  array  with  data  from  temporary  buffer  */ 
vadd_(pix_index,&i_one,pix_temp,&i_one,pix_index,&i_one,&i_pix„size); 

/*  pix . index  jumps  along  k  axis  */ 

pix_index=pix_index+pix_index_inc; 

}  /*  k  loop  */ 


}  /*  subroutine  */ 


B.6.  Declare  AH  Global  Constants  (Focus.h) 

#define  NO_SUPERCARD  3  /*  number  of  supercards  */ 

#define  NPOINT_AZIMUTH  600  /*  number  of  bearing  lines  */ 

#deflne  NPOINT_RADIAL  4095  /*  number  of  radial  lines  */ 

#define  NPOINT_APER  2304  /*  #  of  positions  in  aperture  */ 

#define  NPOINTJDATA  2048  /*  number  of  original  data  points  */ 
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#define  NPOINT_DATA__INTER  NP0INTJDATA*8  /*  number  of  data  points  after 

interpolated  */ 

^define  FILTER  JLENGTH  255  /*  length  of  FIR  filter  */ 

^define  PI  3.141592654 

#defme  RADAR„HEIGHT  60.  /*  radar  height  in  feet  */ 

#deftne  AJ3FSET  - 15.0  /*  offset  angle  from  center  line  */ 

#define  A_SPAN  54.0  /*  spanning  angle  of  the  patch  */ 

#define  RMINCENTER  550.  /*  range  from  center  position  to 

reference  point  on  tha  patch  */ 

#define  RCENTER  REF  802.0  /*  range  (ft)  of  ref.  point  from  center  pos.  */ 

#define  THETA_CENTER_REF  -15.0  /*  angle  (deg)of  ref.  point  from  center  pos.  */ 
#define  SAMPLING_RATE  2.0e09  /*  sampling  frequency  of  signal  */ 

#define  SAMPLING^  PERIOD  5.0e-10  /*  ts=l/fs  */ 


#define  BOX_S  IZE_RADIAL  1 95 
#define  BOX_SIZE_  AZIMUTH  100 

#define  NBOX_RADIAL  NPOINT  RADIAL/BOX__SIZE  ,RADIAL 

#define  NBOX_AZIMUTH  NPOINT_AZMUTH/BOX_SIZE_AZIMUTH 

#define  BOX_BASE_START  2 

#define  NO_KBOX_PROCESS  3 

#defxne  IBOX^START  10 

#define  NOJBOX  PROCESS  1 

#define  SECTION„SIZE  4  /*  #  of  sections  for  one  aperture  */ 

#define  COEF_SIZE  6  /*  #  of  coeffs  for  curve  fit  */ 

#define  MAXDEG  3  /*  maximum  degree  for  poly,  fit  */ 

#define  NPOINT_BOXSAMPLE_RADIAL  5  /*  #  of  samples  in  a  box  in  radial 

direction  */ 

#define  NPOINT  BOXSAMPLE^AZIMUTH  5  I*  #  of  samples  in  azimuth  direction 

for  every  box  */ 

#define  NPOINT_BOXS  AMPLE 

NPOINT_BOXSAMPLE_RADIAL*NPOINT_BOXSAMPLE__AZIMUTH 

/*  #  of  samples  in  a  box  for  curve  fit  */ 

B.7.  Declare  AH  Global  Variables  (Focus.var) 

struct  complex  {  float  r,i; } ;  /*  define  a  complex  type  */ 

struct  cstype  { 

/*  Data  Variables  */ 

float  pixel[NPOINT_AZIMUTH/NO_SUPERCARD][NPOINT_RADIAL];  /*  array  of 
image*/ 

float  data_temp[NPOINT_D ATA] ;  /*  temp  buffer  for  data  */ 

float  data[NPOINT_D ATA] ;  /*  original  data  bufer  */ 


81 


float  filter_coef[NPOINTJDATA_INTER];  /*  filter  coefficients  */ 

float  data_inter[NPOINT_DATA_INTER+136];  /*  interpolated  data  buffer  */ 

float  tstart_err[NPOINT_APER];  /*  tstart-tstart_trunc  for  each  position  */ 

struct  complex  data_inter_fft[NPOINT„DATA„INTER/2];  /*  Real->Complex  Forward  FFT 

of  data inter  */ 

struct  complex  filter„fft[NPOINT_DATA_INTER/2];  /*  Real->Complex  Forward  FFT 

of  filter„coef  */ 

float  nyquist_filter;  /*  value  of  FFT  of  filter  at  Nyquist  point  */ 

float  nyquist_data;  /*  value  of  FFT  of  data  at  Nyquist  Point  */ 

long  ham_flag;  /*  to  indicate  wheter  to  use  hamming  or  not  */ 

float  ham_coef[NPOINT_APER] ;  /*  hamming  weight  coefficients  */ 


/*  Geometry  Variables  */ 

float  a_ofset;  /*  offset  angle  from  the  center  line  */ 

float  a_span;  /*  spanning  angle  of  the  patch  */ 

float  d_theta;  I*  delta  angle  */ 

float  x, theta;  /*  coordinate  of  a  pixel  */ 

float  c_theta,s_theta;  /*  cosine  and  sine  of  theha  */ 


float  d;  /*  distance  from  radar  to  center  position 

positive  to  the  right,  neg.  to  the  left  */ 
float  rcenter,rmincenter,rmaxcenter;  /*  range  from  center  position  */ 
float  rcenter_ref;  /*  range  from  center  pos.  to  ref  point  */ 

float  theta_center_ref;  /*  angle  of  ref.  point  from  center  position  */ 


float  r_ref;  /*  range  from  any  pos.  to  ref.  point  */ 

float  djrcenter;  /*  sampling  distance  from  center  position  */ 

float  dr8;  /*  (sampling  distance)/8  */ 

float  r,rmin,rmax;  /*  range  from  any  position  */ 

float  aperjength;  /*  length  of  aperture  */ 

float  radar_height;  I*  height  of  radar  */ 

float  radar  Jieigbt2;  /*  square  the  height  of  radar  */ 


/*  Radar  Variables  */ 
float  c; 
float  ts; 
float  fs; 


/*  speed  of  wave 
/*  sampling  period  of  signal 
/*  sampling  frequency  of  signal 


*/ 

*/ 


long  n_azimuth; 
long  n_radial; 
long  n_aper; 
long  pix_size; 
long  n_data; 
long  n_data_inter; 

long  n.  data . inter_half; 

long  n_boxsample; 
long  n_coef; 


/*  #  of  points  in  azimuth  direction  */ 
/*  #  of  points  in  radial  direction  */ 

/*  #  of  positions  in  aperture  */ 

/*  #  of  points  in  pixel  array  */ 

/*  number  of  original  data  points  */ 

/*  number  of  interpolated  data  points  */ 
/*  1/2  #  of  interpolated  data  points  */ 
/*  #  of  samples  in  a  box  for  curve  fit  */ 
/*  #  of  coefficients  curve  fit  */ 
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